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abstract       -°«war* 


The  basic  principles  of  the  Galerkin  finite  element 
method  are  discussed  and  applied  to  two  different 
formulations*  one  using  different  basis  functions  and  the 
other  using  the  vorticity — divergence  forrr  of  the  shallow 
water  equations.  lach  f angulation  is  compared  to  the 
primitive  form  of  the  equations  developed  by  lelley  (1976). 
The  testing  involves  a  comparison  of  three  finite  element 
prediction  models  using  variable  size  elements.  Equilateral 
elements  significantly  improve  the  solution  and  are  used  in 
most  of  the  comparisons.  The  formulation  using  different 
basis  functions  produces  poorer  results  than  the  primitive 
formulation.  The  vorticity-divergence  formulation  produces 
superior  results  while  executing  faster  than  the  primitive 
model.  Fcwever,  it  does  require  more  storage  and  the 
relaxation  parameters  are  sensitive  to  the  domain  geometry. 
The  corputer  implementation  for  the  vorticity-divergence 
model  is  discussed  and  the  source  listing  is  included. 
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I.  INTHOrUCTION 

Shuman  fl978)  claims  that  progress  in  numerical  modeling 
cf  the  general  circulation  has  been  to  some  degree  dictated 
in  the  past  "by  the  rate  of  development  in  the  field  of 
computer  technology.  Eowever,  the  limited  ability  to 
Darameterize  the  effects  of  small-scale  processes  in  terms 
of  large  scale  motions  has  been  an  equally  important 
limiting  factor.  Essentially,  the  major  problem  of  numerical 
modeling  of  the  general  circulation  is  simply  that  of 
producing  a  very  long  range  numerical  weather  forecast. 

Certainly  the  equations  used  in  the  models  must  be  more 
sophisticated  to  include  those  physical  processes- which  are 
unimportant  for  a  short  range  forecast,  but  may  become 
crucial  as  the  length  of  the  forecast  is  extended.  Another 
area  where  concentrated  efforts  have  improved  the  forecast 
involves  the  computational  techniques  employed  to 
approximate  and  solve  the  governing  equations  of  the  models. 

The  motivation  behind  this  thesis  is  to  investigate  the 
application  of  a  relatively  new  computational  technique  to 
the  field  of  numerical  weather  prediction.  The  finite 
element  method,  long  established  in  engineering,  has  been 
seriously  considered  only  during  the  past  decade  In 
meteorology.  This  method  has  great  potential  for  application 
in  atmosoheric  Prediction  models. 
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A.   BACKGROUND 

The  rrost  common  numerical  integration  procedure  for 
weather  prediction  has  been  the  finite  difference  method  in 
which  the  derivatives  in  the  differential  equations  of 
motion  are  replaced  by  finite  difference  approximations  at  a 
discrete  set  of  points  in  space  and  time.  The  resulting  set 
of  equations,  with  appropriate  restrictions,  can  then  he 
solved  by  algebraic  methods.  Until  recently,  the  finite 
difference  method  has  been  the  workhorse  in  atmospheric 
prediction  models,  from  their  first  computer  implementation 
to  the  present . 

With  the  introduction  of  each  new  generation  of 
computers.  the  gap  between  numerical  forecasts  and 
atmospheric  observations  has  decreased.  The  rate  at  which 
this  gap  decreased  has  slowed  down  and  appears  to  be 
leveling  off.  This  would  indicate  that  computer  technology 
may  not  be  the  prirrary  obstruction  to  better  numerical 
forecasts.  In  fact,  bigger  and  faster  computers  alone  have 
demonstrated  their  inability  to  significantly  improve  the 
numerical  forecast. 

7cr  example,  a  major  limiting  factor  of  finite 
difference  approximations  is  the  truncation  error.  The 
National  leather  Service  7  Layer  Primitive  Equation  P-odel 
(7LPZ  Model),  operational  from  1966  to  1980.  had  truncation 
errors  which  increased  at  a  rate  proportional  to  the  square 
o*"  the  grid  spacing.  That  is,  the  smaller  the  grid  interval, 
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the  siraller  the  truncation  error.  To  increase  its  accuracy 
would  require  increasing  the  grid  matrix  density.  This  would 
require  increased  computer  storage  and  computational  tiire. 
State  of  the  art  computers  are  capable  of  providing  these 
additional  resources. 

The  problem  now  goes  beyond  numerical  techniques  and 
romputer  technology.  Operationally,  the  National  Weather 
service  is  not  capable  (due  to  monetary  restrictions)  of 
providing  a  denser  concentration  of  atmospheric 
observations.  Therefore,  with  the  present  density  of  initial 
data  ( observations)  and  objective  analysis  techniques 
(getting  the  data  for  grid  points  by  interpolating  from 
observed  data  sources),  reducing  the  grid  spacing  further  on 
the  7LPZ  Model  does  not  significantly  increase  the  accuracy 
of  the  solution. 

This  additional  computer  capability  can  not  be  utilized 
usin«?  finite  difference  methods.  Therefore,  new  numerical 
integration  techniques  must  be  investigated,  such  that  given 
the  same  density  of  observed  data,  superior  solutions  are 
produced . 

Two  alternative  techniques,  the  spectral  method  and  the 
finite  element  method,  have  started  to  gain  attention.  Both 
the  spectral  and  finite  element  methods  require  more 
computational  time  per  forecast  time  step  than  does  the 
finite  difference  method.  For  example,  the  finite  element 
method  requires  an  equation  solver  to  invert  a  larger  matrix 
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at  each  tirre  step  for  each  variable.  In  this  sense,  these 
methods  were  held  hack  by  computer  technology,  hut  recent 
advances  in  computer  technology  (i.e.  larger  and  faster 
storage  devices,  multi-processors,  etc)  have  made  these 
alternative  numerical  techniques  competitive. 

For  long  range  weather  predictions,  the  spectral  method 
applied  over  the  globe  or  hemisphere  is  a  natural  method, 
due  to  the  existence  of  efficient  transforms  for  the 
nonlinear  terms  on  spherical  geometry.  It  also  eliminates 
the  truncation  error  for  the  horizontal  space  derivatives 
and  the  nonlinear  instability  (aliasing).  For  these  reasons, 
global  spectral  models  have  been  developed,  and  implemented 
on  an  operational  level,  replacing  the  global  finite 
difference  models . 

Eowever,  because  the  spectral  harmonics  are  globally 
rather  than  locally  defined,  it  is  thought  that  for  problems 
of  more  detailed  limited  area  forecasting,  the  finite 
element  method  is  more  suitable.  Pioneering  work  to  adapt 
finite  element  methods  to  meteorological  applications  has 
been  done  by  Cullen  (1973,1974  and  1979),  Staniforth  and 
Mitchell  '1977),  Kinsman  (1975)  and  lelley  (1976).  The  most 
recent  finite  element  meteorological  model  at  the  Naval 
Postgraduate  School  was  written  by  Kelley  (1976)  with  the 
-?llaboration  of  Tr.  R.T.  Williams.  It  is  this  study  that 
will  serve  as  a  basis  for  this  thesis.  The  model  written  by 
Kelley  will  be  altered  and  used  for  comparative  testing  with 
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irproved  finite  element  fortrs  implemented  by  this  author. 
Some  of  the  techniques  and  coles  developed  by  Keliey  are 
also  employed  in  this  thesis.  Older  (1981)  developed  a 
technique  to  smoothly  vary  the  grid  geometry  in  the  domain. 
This  technique  is  also  implemented  both  on  Kelley's  model 
and  vith  the  new  formulation  to  give  greater  versatility 
when  testing  the  model  performance. 

3.   C3J2CTI7ES 

The  objectives  for  this  thesis  can  be  divided  into  two 
categories:  1)  meteorology,  2)  computer  science.  First,  the 
meteorological  objectives  of  developing  improved  finite 
element  forms  for  shallow  water  equations  are  as  follows: 

1)  -  Older  (1981)  after  collaboration  with  Dr. 
M.J. P.  Cullen,  showed  how  equilaterally  shaped  elements 
produced  significantly  better  results  than  did  other 
triangular  elements,  ^elley  (1976)  used  right  triangular 
elements  in  the  implementation  of  a  two-dimension.al  finite 
element  model  using  the  primitive  form  of  the  shallow  water 
equations.  A  considerable  amount  of  small-scale  noise  was 
observed  in  the  solution.  Hereafter,  this  model,  which  was 
developed  by  lelley  (1975),  will  be  referred  to  as  the 
primitive  model.  This  first  objective  involves 
re-implementing  the  primitive  model  using  equilaterally 
shaped  elements  and  comparing  the  results  to  those  in 
Zelley's  thesis. 
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2^  -  villiaT>s  and  Zienlciewicz  (1981)  presented  new 
finite  elerrent  techniques  for  formulations  for  the  shallow 
water  equations,  which  use  differently  shaped  functions  to 
approximate  the  different  dependent  variables,  which  in 
effect  stagger  the  variables.  Schoenstadt  (1980) 
demonstrated  the  advantage  of  spatial  staggering  of 
dependent  variables  in  finite  difference  models.  The 
application  of  this  technique  to  finite  element  models  is  a 
natural  extension,  and  excellent  results  were  obtained  by 
Williams  and  Zienkiewicz  (1981)  from  application  of  these 
r.ew  formulations  on  linearized  one  dimensional  cases.  The 
objective  here  is  to  implement  the  new  forms  on  the 
primitive  jodel  and  again  do  quanitative  comparisons  of  the 
results . 

?)  -  The  major  emphasis  in  this  study  deals  with  the 
implementation  and  comparison  of  the  vorticity  divergence 
for'r  of  the  shallow  water  equations,  which  is  described  in 
detail  in  Chapter  III.  This  formulation  has  the  following 
advantages.  First,  the  geostrophic  adjustment  process  is 
treated  better  than  in  the  primitive  form  of  the  equations. 
Secondly,  the  velocity  and  height  fields  are  evaluated  at 
the  same  grid  point,  where  the  best  primitive  form  requires 
staggering  these  dependent  variables.  And  thirdly,  a  larger 
time  step  is  allowed  due  to  the  semi  implicit  form  of  the 
equation.  Again  comparisons  between  the  results  from  the 
vorticity  divergence  and  primitive  model  are  presented. 
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The  computer  science  aspect  of  this  thesis  was  primarily 
devoted  to  the  implementation  of  the  different  models  and 
the  style  and  architecture  of  the  program.  Finite  element 
methods  require  more  computational  time  than  do  finite 
difference  methods,  not  only  in  the  solution  of  the 
equations,  but  also  in  the  amount  of  computation  required  to 
evaluate  each  term  in  the  equations. 

The  implementations  of  these  two  dimensional  models, 
although  complex  when  viewed  from  the  surface,  have  a  lot  of 
generality  and  redundancy  in  the  operations  required. 
Versatile  modules  can  he  written  to  ease  the  implementation 
and  facilitate  changes.  The  objective  here  is  to  efficiently 
implement  these  new  forms  and  demonstrate  tne  utility  of 
these  versatile  modules  for  future  implementations. 

C.   THZSIS  STRUCTURE 

This  thesis  presents  the  results  obtained  from  tests  of 
the  various  finite  element  formulations.  The  results  are 
compared  to  those  from  the  primitive  model  written  by  Celley 
(1976).  Accompanying  the  results  is  a  detailed  discussion  of 
the  ref ormMiation  and  implementation  process. 

Chapter  II  of  the  thesis  presents  a  tutorial  of  the 
finite  element  method  and  the  area  coordinates  system  used 
in  the  evaluation  of  the  element  integration.  The  Galerkin 
finite  element  method  used  in  all  the  models  is  developed 
and  applied  to  the  advection  equation  in  one  dimension. 
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rhapter  in  presents  the  detailed  description  cf  the 
vcrf i ri ty-di vergence  shallow  water  model.  Here  the  equations 
are  shown  and  written  using  the  Salerkin  method.  A 
discussion  of  the  computational  technique  used  is  presented 
along  with  the  -model's  physical  parameters. 

Chapter  IV  presents  a  descriptive  overview  of  the 
computer  implementation.  The  chapter  includes  a  list  of 
options  available  for  testing,  a  brief  description  of  the 
iratrij  compaction  technique  and  the  formulations  of  the 
versatile  modules  used  to  implement  the  complex  equations. 

Chapters  V  through  VII  discuss  the  results  obtained  fror 
the  different  experiments.  Chapter  V  briefly  describes  the 
primitive  model  used  for  all  comparisons  and  the  results 
froffl  changing  the  element  shape  to  equilateral  triangles. 
Chapter  VI  reformulates  the  primitive  model  so  that  the 
geopotential  is  staggered  with  respect  to  the  velocity 
variable.  Jor  simplicity,  the  continuity  equation  is  also 
linearized.  Chapter  VII  compares  the  results  from  the 
vert icity-di vergence  model  developed  in  Chapter  III  to  those 
:"rom  the  primitive  model. 

The  last  chapter  summarizes  the  results  from  all  the 
experiments  and  identifies  what  areas  require  follow  on 
work.  The  source  code  for  the  vorticity-divergence  model  is 
presented  in  Appendix  A. 
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II  -  FIN  ITS  5LZMINTS 

As  is  often  the  case  with  an  original  development,  it  is 
rather  difficult  to  quote  an  exact  date  on  which  the  finite 
element  method  was  invented,  but  the  roots  of  the  method  can 
be  traced  tack  to  these  separate  groups:  applied 
mathematicians,  physicists  and  engineers.  Since  the  early 
developments  of  the  finite  element  method,  a  large  amount  of 
research  has  been  devoted  to  the  technique.  However,  the 
finite  element  method  obtained  its  real  impetus  by  the 
independent  developments  carried  out  by  engineers.  Its 
essential  simplicity  in  both  physical  interpretation  and 
mathematical  form  has  undoubtedly  been  as  much-  behind  its 
popularity  as  is  the  digital  computer  which  today  permits  a 
realistic  solution  of  even  the  most  comolex  situations. 

The  name  finite  element  was  coined  in  a  paper  by 
R.V.  Clough,  in  which  the  technique  was  presented  for  plane 
stress  analysis,  as  discussed  in  Bathe  (1976).  While  finite 
element  methods  have  made  a  deep  impact  via  the  field  of 
solid  mechanics,  where  it  can  be  said  that  today  they 
represent  the  generally  accepted  method  of  discretizing 
continuum  problems  for  computer-based  solution,  the  same 
appears  not  to  be  true  in  fluid  mechanics  or  atmospheric 
ored  ic tion  . 


19 


Numerous  finite  element  formulations  are  currently 
available.  Strang  (1973),  Ncrrie  (1973)  and  Zienkiewicz 
11571}  present  detailed  theoretical  discussions  of  each.  The 
Galerkir.  method  ,  the  most  popular  finite  element  method,  is 
described  in  detail  below  and  used  in  the  equation 
formulation  later. 

A.   FINITE  ILZMINT  CONCEPT 

The  problem  of  solving  partial  differential  equations 
can  be  specified  in  one  of  two  ways.  In  the  first,  finite 
difference  methods  specify  the  dependent  variables  at 
certain  grid  points  in  space  and  time,  and  the  derivatives 
are  evaluated  using  Taylor  series  approximations.  Secondly, 
the  calculus  of  variation  requires  the  minimization  of  a 
functional  over  a  domain,  where  a  functional  is  defined  as  a 
variational  integral  over  the  domain. 

The  calculus  of  variation  approach  creates  a  purely 
physical  model  where  the  functional  equivalent  to  the  known 
differential  equations  are  lenown.  Its  major  disadvantage  is 
that  it  limits  the  method  only  to  those  problems  for  which 
fun'Uionals  exist.  Finite  element  methods,  an  extention  of 
this  method.  derive  mathematical  approximat icns  directly 
fror  the  differential  equations  governing  the  problem.  The 
advantage  here  is  that  it  extends  the  method  to  a  range  of 
problems  for  which  a  functional  may  not  exist,  or  has  not 
been  ii  sc  ove  red  . 
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The  finite  eleirent  rrethod  divides  the  domain  into 
subdorains  or  finite  elements  (usually  of  the  sarre  form). 
v.Mes  are  located  along  the  boundary  of  the  elements, 
usually  at  the  element  vertices  and  at  strategic  positions 
(mid! side,  centroid,  etc.)  in  the  interior  and  on  the  sides 
of  faces  of  eletrents. 

Commonly  used  elements  are  triangular,  polygonal  or 
polyhedral  in  form  for  two-dimensional  problems.  The  choice 
of  elements  depends  on  the  type  of  problem,  the  number  of 
elements  desired,  the  accuracy  required  and  the  available 
computing  time.  To  begin  with,  the  element  must  be  able  to 
represent  derivatives  of  up  to  the  order  required  in  the 
solution  procedure,  and  to  guarantee  continuous  first 
derivatives  across  the  element  boundaries  to  avoid 
singularities.  Triangular  elements  are  errployed  in  this 
thesis  because  they  can  be  used  effectively  to  represent 
irregular  boundaries,  and/cr  geometry,  and  also  to 
concentrate  coordinate  functions  in  those  regions  of  the 
domain  where  rapidly  varying  solutions  are  anticipated. 

Consider  the  problem  of  solving  approximately  the 
differential  equation 


L(u)  =  f(x 


II-l 


where  I  is  a  differential  operator,  u  the  dependent 
variable,  and  f(x)  is  a  specified  forcing  function.  Suppose 
that   II-l   is  to  be  solved  in  the  dorrain  a  *  x  -  b  and  that 
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appropriate  boundary  conditions  are  provided.  The  residual  P. 
is  forced  fror  II— 1  as  follows: 


L(u)  -  f(x)  =  R 


II-2 


The  critical  step  is  to  select  a  trial  fairily  cf 
approximate  solutions  (the  members  of  a  trial  family  are 
often  called  basis  functions).  The  basis  function  is 
prescribed  'functionally)  over  the  domain  in  a  piecewise 
fashion,  element  by  element,  and  are  generally  a  combination 
of  low  order  polynomials .  A  one  dimensional  example  is 
shown  in  Figure  1,  wherein  the  domain  (x  axis)  is  divided 
into  six  elements  (line  segments)  A  through  F.  The  basis 
functions  are  linear  and  one  is  shown  for  node  4  only  in 
Figure  1.  mhe  function  has  a  value  of  unity  over  node  4,  and 
decreases  linearly  to  is  zero  at  nodes  3  and  5  and  zero 
els  swhere. 

Consider  a  series  of  linearly  independent  basis 
functions  V.  (x),  as  in  Figure  1.  Now  u(x)  can  be 
approximated  with  a  finite  series  as  follows: 


u(x)  =   Z  $.  V. (x)  =  0.V. 

\     J  J       J  J 


where  £  is  the  coefficient  of  the  jth  basis  function  and  has 
j 

a  value  equal  to  u  at  node  j. 

Substituting   this   approximate  solution  I 1-3  wherever  u 
appears  in  the  differential  equation  1 1  — 1 


L^.V.)  -  f(x)  =  R 
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The  best  solution  will  be  one  which  in  some  sense 
reduces  the  residual  R  to  a  minimum  value  at  all  points  in 
the  domain.  By  definition,  the  residual  obtained  using  the 
exact  differential  equation  is  identically  zero  everywhere. 
The  residual  P.,  formed  in  equation  II-4,  is  minimizied  when 
multiplied  with  a  weighting  function,  integrated  over  the 
domain  and  set  equal  to  zero.  This  process  is  known  as  the 
weighted  residual  method 


b 

RV  dx  =  3  II-5 


\ 


where  W  is  the  weighting  function  and  is  referred  to  as  the 
test  function  in  the  following  development.  The  weighted 
residual  method  minimizes  the  errors  of  the  residuals,  such 
that  the  summation  of  all  the  positive  and  negative  errors 
add  to  zero . 

The  ^alericin  method,  the  most  popular  fini  te  element 
method,  is  more  general  in  application  and  is  a  special  case 
of  the  method  of  weighted  residuals,  as  discussed  by  Pinder 
ar.d  Gray  ^1977).  The  requirement  imposed  on  the  weighted 
residual  method  forming  the  ^alerkin  method  is: 

*  the  test  (weighting)  function  be  equal  to  the 
basis  (trial)  function  W  =  V  .  This  process 
leads  in  general  to  the  best  approximation  of 
the  solution. 
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The   final   (Jalerkin   form   is   obtained   by   substituting   II-4 
into    II-5,    yielding 


-fi±l((j}.l.  )dx  -         Vf1f  (x)dx   = 


II-S 


If  this  procedure  is  repeated  for  N  points  in  the  domain 
a  system  with  N  equations  and  N  unknowns  will  be  generated. 

3.   OAimiN  APPLICATION 

The  following  example  taken  in  part  from  Haltiner  and 
Williams  (1980)  applies  the  Galerkin  rethod  to  the  advection 
equation  with  linear  elements 


du     du 

—  +  c  —  =2 

at        dx 


II-7 


This  equation  is  dependent  in  both  time  and  space.  The 
treatment  of  time  variation  is  important  for  most 
meteorological  prediction  problems.  The  C-alerkin  method  is 
not  applied  to  the  time  dependence  because  it  is  more 
convenient  to  use  finite  differences  in  timefas  is  done  with 
this  example  later.  The  same  treatment  is  applied  tc  the 
prognostic  equations  later,  where  two  finite  differencing 
methods  are  employed  to  do  the  time  integration. 

The  7alerkin  procedure  represents  the  dependent  variable 
u'x.t)   with   a   sum   of   functions  that  have  the  prescribed 
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spatial  structure  as  in  Figure  1.   Approximate   u(x,t)   with 
the  finite  series  as  follows 


u(xft)  =  I   rf.(t)V .(i)  =  0.?, 

1  =  1^    ^       ^  *J 


II-8 


where  the  coefficient  0  .  (t),  a  function  of  tirre,  is  the 
scalar  value  of  u  at  node  j.  The  basis  functions,  V.(x),  are 
functions  of  space  only  and  j  equals  1  to  ?  for  the  example 
In  Figure  1.  The  repeated  subscript  in  this  forx  implies  a 
sun-  ever  the  repeated  subscript. 

The  Galerkin  equation  for  the  advection  equation  II-7  is 
obtained  by  setting  L  =  c(d()/dl)  and  substituting  in  the 
aooroximate  solution  II- 8  wherever  u  is  found. 


N 
1 

j=i  at 
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m    J 
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— JV, 
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II-9 


where  i  =  1  to  Nv  V.  the  test  function  and  V.the  basis 
function.  The  domain  of  integration  is  given  by  a  *  x  -  b, 
and  the  integration  is  done  in  a  piecewise  fashion,  element 
by  element . 

In   this  one-dimensional  case,  an  equation  like  1 1—9  is 
written   for  each  node,  i.  Considering  node  4,  what  are   the 
possible  non-zero   contributions  from  equation  1 1—9?  Figure 
2  illustrates   the  basis  and  the  test  function   interaction 
luring   the    piecewise   integration   process.   From   the 
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Figure  2 .   Basis  and  test  function  interaction 
during  the  piecewise  integration  process . 
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definition  of  the   basis  and   the   test   function,   locally 
defined   as  unity  at  node  j  and  linearly  decreasing  to  zero 

;  ±  1  and  zero  elsewhere,  the  only  non — zero 
contributions  are  made  when  *  =  3  over  element  C,  j  -  4  ever 
elements  ~  and  ?,  and  j  =  5  over  element  D. 

The  evaluation  of  1 1-9  for  i  =  m,  which  is  given  in 
Haltinep  and  Williams  (1981),  leads  to  the  equation: 


6  dt    m+1 


4u   +  u  1  )  +   e  ( u   1  -  um  -  )  =  0    11-10 
in     m-1     0     m+1     m-1 
2-*x 


The  boundary  points,  which  in  this  example  are  nodes  1 
and  7,  are  evaluated  in  the  same  way  as  the  interior  nodes, 
with  the  exception  that  cyclic  conditions  are  imposed. 

The  time  discretization  of  II— 10  is  d-one  using-  a  finite 
difference  scheme.  Applying  leapfrog  time  differencing  gives 
the  following  equation 


f 


12  At 


m+1     m+1       m      m 


n+1 
m-1 


n-1  , 
m-l 


Hi    m+1 


)  =  0 


n-ii 


The  resultant  equation  set,  in  matrix  form,  contains  an 
NxN  matrix  where  N  is  the  number  of  nodes. 

The  transition  from  one-dimension  to  two  is 
mathematically  identical.  The  domain  is  now  subdivided  into 
finite  areas,  which  are  triangles  in  this  implementation  and 
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the  basis  functions  are  linear.  However,  now  they  are 
pyrarid  shaped  with  value  unity  at  the  center  and  decrease 
to  zero  at  the  surrounding  nodes,  and  are  zero  elsewhere. 
Figure  2  shows  this  basis  function  for  node  28  outlined  in 
heavy  black.  The  value  at  any  node  again  can  be  approximated 
by  II-3,  where  j  ranges  over  all  nodes  connected  to  node  i 
including  i  itself.  The  connectivity  for  node  i  =  28  in 
Figure  2  is  j  =  15,16,27,28,29,39  and  40. 

The  integration  is  still  over  the  entire  domain.  With 
both  the  basis  and  the  test  function  zero  over  the  domain, 
except  locally  over  each  element,  the  global  integration  can 
be  performed  by  integrating  locally  over  each  element.  By 
definition,  this  integration  can  be  expressed  as  an  inner 
product  of  both  functions  (i.e.  basis,  test)  as  follows: 


<w  = 


re 
)) 

A 


V .  V,  dA 


Using   this  definition  and   the   repeated 
notation  eauation  II— 9  becomes 


11-12 


subscript 


»-<?„?!>  ♦  c  iM*.X<>   =  0 


11-13 


where  the  dot  implies  differentiation  with  respect  to  time, 
and  the  second  subscript  implies  differentiation  with 
respect  to  the  second  subscript.  The  local  integration  may 
be  calculated  directly  from  exact  expressions  derived  from 
area  coordinates  described  in  detail  in  the  next  subsection. 


C6 


Figure  3.  Basis  function  for  node  28.  The 
shaded  area  is  the  complete  basis  function 
and  the  V.  ,  where  j  =  15,  16,27,28,29,39,^0 
are  jth  node  basis  functions  for  node  28.  The 
dashed  line  at  node  28  has  length  unity. 
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In  s unwary,  the  Galerkin  procedure  involves  subdividing 
the  domain  into  finite  elements,  approximating  the  dependent 
variables  by  a  linear  combination  of  low  order  polynomials 
a~A  substituting  them  into  the  equations.  The  equation  is 
multiplied  by  a  test  function,  integrated  over  the  entire 
domain  and  finally  the  resulting  system  of  equations  is 
solved . 

C.   APIA  COORDINATES 

While  the  Cartesian  coordinate  system  is  the  natural 
choice  of  coordinates  for  most  two  dimensional  problems,  it 
is  not  convenient  when  working  with  triangularly  shaped 
elements.  It  is  therefore  necessary  to  define  a  special  set 
of  normalized  coordinates  for  a  triangle.  Area,  or  natural 
coordinates  as  they  are  commonly  called,  reduce  the 
formidable  task  of  integrating  products  between  the  basis 
avA  test  functions  and  their  derivatives  over  a  triangular 
element  and  result  in  easily  computable  and  exact 
expressions . 

The  following  development  is  taken  in  part  from  the 
formulation  by  Zienkiewics  (1971).  Consider  the  triangular 
element  illustrated  in  Figure  4.  There  is  a  one-to-one 
correspondence  between  the  Cartesian  coordinates  (X,Y)  and 
the  area  coordinates  (L^Lg.L-a  )  for  the  element.  Let  A 
denote  the  area  of  the  triangle  and  Aj_  ,  ^and  A^  the  areas  of 
the  svbtriangles  in  Figure  4  such  that  A  =  A^+Ag+A,. 
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Fig.  5  •  Transformation  to  area  coordinates 
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The  relationship  between  a  point  P(X,I)  in  Cartesian 
rcordinates  and  ?(L  ,  L  ,  L_  )  in  area  coordinates  can  be 
seen  by  the  following  transformations 


X  =  LjX   +  L2X   *  L-X 

Y  =  L1Y   +  L^Y   +  L3Y 

1  =  Ll  +  L2  +  L3 


11-14 


where 


■.-*■ 


L-  =  -  and 
*   A 


and 


L1  =  f2A  +  b1X  +  a1Y)  /2A 
L2  =  (2A  +  b2X  +  a2Y)  /2A 
L3  =   '2A  +  b3X  +  a3Y)  /2A 


11-15 


where  2A   is  twice  the  area  of  the  triangle  and  the  a's  and 
b's  are  defined  as  in  Figure  5. 

It  is  worth  noting  that  every  tuple  (L-  ,  L2  ,  L-»  ) 
corresponds  to  a  unique  pair  fX,Y)  of  Cartesian  coordinates. 
In  Figure  4,  L-a  1  at  vertex  1  and  0  at  vertices  2  and  3.  A 
linear  relation  exists  between  the  area  and  Cartesian 
coordinates  which  iirplies  that  values  for  L.vary  linearly 
over  the  triangle  with  a  value  one  at  vertex  1  and  a  value 
of  zero  at  vertices  2  and  3;  and  similarly  for  L2and  L- . 
This  demonstrates  how  each  component  in  the  tuple  (L^,  L2 ,  LJ 

behaves  over  the  triangle  as  do  the  linear  basis  and      test 
functions  over  the  element,  as  was  seen  in  Figure  4.  Clearly 


=  7. 


11-16 


33 


where   V.  is  a  linear  function  of  the  Cartesian  coordinates 
(i.e.  "basis  ,  test )  . 

Zienkievicz  (1971)  shows  that  it  is  possible  to 
integrate  any  polynomial  in  area  coordinates  using  the 
sirole  relationship 


~  1  ^2  ^3  dxd y  = 


m !  n  !  o  ! 


A 


( m.  +  n  +  p  +  2 ) ! 


2A 
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where  rrf  a  and  p  are  positive  integers  and  A  is  the 
elementary  area.  For  an  example  of  this  integration 
technia.ue  using  inner  product  notation,  equation  11-12  is 
evaluated  as  follows 


V^  didy  = 


2!  0!  0! 


2A  =  - 

(2  ♦  0  +  0  +  2)!       6 


<vv-< 


f 


1!  1!  0! 


V.V,  dxdy  =  

J  (1  +  1  +  0  + 


A 

2A  =  — 

2)!       12 


i  =  J 


i  t   J 
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The  differential   operations  in  area  coordinates  follow 
directly  from  the  differentiation  of  (11-15)  where 


ar.t 


b     3  b.   o 

bx   i=l  2A  dL4 


b    3  *      a 

'by        i=l  2A  bL. 


11-19 


11-20 
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As  explained  earlier  (see  Equation  11-16),  Viis  a  linear 
function  'i.e.  basis,  test)  which  equals  a  component  L1  of 
the  srea  coordinate  tuole.  Therefore 


5  = 


0  if   1  ^  j 

1  if   i  =  j 
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Consequently   BV.  dx 


jx 


dv. 
dx 


h.  I!1 

2  A   dl< 


for  j  =  1   is 


2A 
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Is  an   exarrple,   consider  the  inner  product  <V  .  ,v.>  at 
vertices  j  =  2,  i  =  1.  This  integration  is  evaluated  as 


<72x»V  = 


(f  b 


y 


>  2A 


71   ±xdy 
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.!* 


1!  0!  0 


2A   =  - 

2A   (1+0+0+2)!        6 


Therefore  any  inner  product  in  the  formulation  can  be 
readily  evaluated  using  area  coordinates.  Another  benefit  of 
using  this  coordinate  system  is  that  all  of  the  inner 
products  are  functions  of  space  only  and  need  be  computed 
only  once. 
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III.   SHALLOW  WATIR  ^Clil 


The  governing  equations  for  this  model  are  derived  by 
raking  several  simplifying  assumptions  on  the  primitive 
equations  of  motion,  whish  then  give  the  barctropic  shallow 
water  equations.  However,  as  mentioned  previously,  the 
shallow  water  equations  describe  many  significant  features 
of  the  large-scale  motion  of  the  atmosphere,  and  therefore 
have  been  used  in  numerous  experiments  over  the  years. 

The  vcrticity-divergence  form  of  the  equations  has 
several  advantages.  Williams  (1981)  has  shown  that  the 
frecstrophic  adjustment  process  is  treated  much  better  with 
the  vcrticity  divergence  formulation  than  with  a  direct 
treatment  of  the  primitive  form  of  the  shallow  water 
equations,  such  as  was  used  by  Kelley  (1976).  This 
formulation  also  allows  the~  velocity  components  and  the 
height  to  be  carried  at  the  same  nodal  points,  whereas  the 
best  scheme  for  the  primitive  form  of  the  equations  requires 
staggering  of  the  fields,  as  seen  in  Schoenstadt  (1980).  The 
▼orticity  divergence  form  of  the  equations  is  also 
convenient  for  the  application  of  semi-implicit 
differencing,  which  saves  considerable  computer  time. 


A.   SOYIRNING  EQUATIONS 

The   primitive   form   of   the  shallow  water  equations  in 
Tirtesiac  coordinates  is 


dt  dx  dy 

du  30  dK 

dt  ox  dx 

dt  dy  dy 


III-l 


III-2 


III-3 


equation  (III-l)  is  the  continuity  equation   and   the   III— 2 

ani      1 1 1-3   are   the  momentum  equations,  respectively.  The 
variables  are  defined  as  follows: 
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the  spatial  coordinates  o*  the  domain 
components  of  the  wind  vector 

geopotential  =  (gravity  x  free  surface  height) 

2         2 

mean  geopotential   =  49,000  meters  /seconds 

time 

kinetic  energy 

absolute  vorticity  =  (-9  +  fj 

relative  vorticity 

ccriolis  force  (mid-channel  f -plane) 

divergence 


The   shallow  water   equations    can 
vorticity  divergence  form  as  follows: 
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where  III -4  is  the  sare  continuity  equation  as  III— 1 ,  I I 1—5 
is  the  vorticity  equation  and  III— 6  is  the  divergence 
equation . 

3ecause  of  the  vorticity  divergence  fortr  of  the 
equations,  it  becomes  necessary  to  solve  the  tirre  dependent 
variables  "S  and  I  in  terms  of  4/,  the  strearr  function 
(rotational  part  of  the  wind),  and*X,  the  velocity  potential 
(divergent  part  of  the  wind).  The  initial  fields  for  the 
model  will  be  in  terms  of  ¥,  TL   and  $, 

The  following  diagnostic  relationships  are  defined  and 
used  later  in  the  solution  of  the  equation  set. 


u  =  "V*x> 


III-7 
III-8 


where  the  subscript  implies  differentiation, 

2    2 

u  ■*■  V 

K   =     kinetic   energy,  I II— 9 

uC      =     u(*    ♦  fj,  111-10 

vC      =      v(S    +  fj  ,  III-ll 

Of  =      0u,  III    12 

P   =      pv,  I I 1-13 
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-8  -  v2  V,  111-14 


D  -   v2*.  111-15 


3.   IOUATION  FORMULATION 

The  G-alerkin  method  described  in  Chapter  II  is  now 
applied  to  eauations  III-4  through  111-15.  For  ease  of 
comprehension »  the  shorthand  inner  product  notation  as  in 
11-12  will  be  used  to  simplify  the  equations.  The  detailed 
Galerkir.  formulation  will  be  shown  for  equation  III-?,  the  u 
component  of  motion.  The  method  follows  directly  from  the 
example  in  Chapter  II  of  this  thesis,  which  in  turn  follows 
in  part  from  £elley  (1976)  and  Ealtiner  and  Williams  (1981). 

Consider  equation  I I 1-7  and  assume  that  each  variable  u, 
^  and  1   is  approximated  by 

-  -   W 

*  =  ¥.  v.  ,  in-ie 

where  the  repeated  subscripts  indicate  summation  over  the 
range  of  the  subscript.  Substituting  these  approximate 
solutions  into  1 1 1-7  yields 

»T   =*   ~^J.)  +  — (*  7.)  1 1 1-17 

Since  only  the  basis  function  7.  is  a  function  of  space, 
TII-17  may  be  further  simplified  by  factoring  out  the  time 
dependent  coefficients. 


The  next  step  requires  multiplying  by  a  test  function   V. 
as   discussed   in  Chapter  II,  and  integrating  over  the  area 
domain 


(   av. 

— JV,  dA 

ay    i 


v.  v,  dA  =  -  ¥ 


A 


J  J 


((    5V. 

— JV,  dA 
dx   x 


III-1S 


'he  final  form  in  inner  product  notation  is 


<uj  Tj-Ti>  = 


"    «J7JT'T1>     +    <V±C'Ti>  nI-1S 


J      JX 


where  the   double   subscript   implies   differentiating  with 
respect  to  the  second  subscript. 

The   three  prognostic  equations  (III-4,  III— 5  and  III-6) 

are  similarly   advanced   using   the   Galerkin   technique   to 
become,  respectively: 

<VrTi>  +  l  <rjW  =  -  <^vjx«vi>  '  <*jTjrV     IXI"20 
<tjfj.fi>  =  -<(ug).Vvi>-<(vc).Vvi>in-2i 

,2, 


<r.V.f7,  >  -  <0.V*V  .,7,>  =  <(vQ).V.  ,V,>  -  <(uC)  .V  .  ,T,  > 
J  J   i      *J    J   i  J  Jx   i  J  jy   i 


+  <£.v27  .,V.  > 
J    J   i 


TT  T-?? 


where  v  is  the  Laplacian  operator  and  the  dot  implies 
differentiation  with  respect  to  the  time  dependence  in 
III-4,  III-  5  and  III-6. 

Similarly,    Galerkin    equations    are   formulated   for 
equations  III-7  through  1 11-15 . 
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C.   TIMZ  riSCP.ETIZATION 

The  equation  set  1 1 1-20,  111-21  and  111-22  is  arranged 
so  that  all  the  tern's  on  the  left  hand  side  can  "be  treated 
Implicitly,  and  all  the  terms  on  the  right  hand  side  can  he 
treated  explicitly.  The  explicit  time  integration  will  be 
done  by  the  leapfrog  difference  method.  To  start  the  time 
integration,  two  forward  half  steps  are  taken,  after  which 
the  full  leapfrog  scheme  is  used  for  the  remainder  of  the 
forecast  period. 

The  vorticity  equation  III- 21  is  solved  independently 
frcr  111-20  or  111-22.  However,  111-20  and  111-22 
'continuity  and  divergence  equations,  respectively)  are 
coupled.  To  explicitly  solve  either,  decoupling  of  the 
equations  is  necessary.  In  this  thesis  this  is  done  through 
algebraic  substitution  of  II 1-22  (solved  for  D(n+1))  into 
IIT-20.  Once  the  time  integration  is  performed  on  111-20, 
111-22  can  he  solved  'or  E(n+1)  using  the  i    (n+1)  value. 

The  final  prediction  equations  are 


-  [3EHT]n+1  -  [BEEYJ  °~1 

"  2Ej[<Tjx-Vix>  +  <VjyViy>] 


111-23 


whe-e   A   =   4/(2*t)    f    E    =  A/J    f    C   =   B/(2*t)    and    [BERY]    is      the 

geostrcphic    boundary   contribution,    see   Section   3. 


.n+1 


v?j.f1> 


n-1 
*j      <Vj,Vi> 


J3 


(^Wi>] 


I I 1-24 


rj    T'j'7i'  - 


DnJ1<7jf7i>    -    (At/2)[0nj+1<v2vj,Vi> 
*   ^S'V    *   2(TQ)3<fja,?1> 


111-25 


After  these  three  elliptic  equations  are  solved,  the 
history  of  the  variables  I I 1-7  through  111-15  is  updated. 

A  large  time  step  can  he  applied  to  this  forrr  of  the 
shallow  water  equations  due  to  the  semi-implicit  nature  of 
the  equations,  -his  is  very  important  since  finite  element 
methods  generally  require  more  computer  time  per  time  step. 
The  vcrticity-divergence  formulation  acts  as  a  filter,  which 
slows  down  the  high  frequency  waves  in  the  solution.  The 
two-dimensional  advective  stability  criterion  for  a  linear 
element,  derived  by  Cullen  (1973),  was  used  to  determine  the 
correct  time  step, 

4  x 


^t  = 


Id  i/e 


111-26 


whereat  is  the  time  step  in  seconds,  ax   the   shortest   grid 
spacing  in  meters  and  c  the  fastest  phase  velocity. 


42 


D.   COMPUTATIONAL  TECHNIQUES 

The  final  prognostic  equation  set  requires  the  solution 
of  a  Felmholtz  equation  for  $  and  Poisson  equations  for  ¥ 
and  J..  The  rrost  common  method  of  solution  used  by 
meteorologists  has  been  the  successive  over  relaxation 
method  (SOH)  in  which  an  initial  guess  of  the  solution  is 
trade  and  then  progressively  improved  until  an  acceptable 
level  of  accuracy  is  reached.  SOH  is  employed  in  the 
solution  of  the  equations, where  111-23  can  be  represented  by 

v2[M]{x}  -  C[M]{x}  *   {b}  111-27 

and  III  24,  III  25  by 

v2[M]{x}   =   {b}  1 1 1-28 

2 

where  v    the  Laplacian  operator,  [M]  =  <V-,V.>  matrix,   {x} 

-  the  dependent  variable  in  vector  notation,  C  -  constant  as 
in  III -23  and  {b}  the  right  hand  side  of  the  equation  or  the 
forcing  function. 

The  mass  matrix  [M] ,  dimensioned  (nxn),  is  a  matrix  of 
coefficients  whose  rows  are  the  equations  of  the  system  to 
be  solved.  There  exists  a  one  to  one  correspondence  between 
the  rows  of  the  mass  matrix  and  the  nodes  of  the  domain. 
Each  equation  has  a  term  (column)  for  each  node,  where  a 
non-zero  term  represents  connectivity.  Nodes  are  connected 
if  they  are  both  vertices  of  the  same  element.  Obviously  [M] 
is   a   sparse   matrix   containing  the  inner  products  for  the 
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left  hand  side.  Chapter  4  of  this  thesis  will  describe  the 
matrix  compaction  procedure. 

The  forcing  function  {b},  dimensioned  (nil),  involves 
only  variables  at  the  current  time  step  and  is  easily 
computed  using  four  very  versatile  subroutines  described  in 
detail  in  the  next  chapter. 

The  initial  guess  to  start  SOR  is  the  previous  time  step 
solution.  An  average  of  30  passes  per  equation  are  needed 
for  each  time  step.  The  solution  is  considered  to  have 
converged  to  its  final  value  when  the  residual  for  each  node 
has  beer,  reduced  to  some  acceptably  small  value. 

The  diagnostic  equations  III-7  through  111-15  must  also 
be  solved  every  time  step.  Eowever,  the  same  technique  is 
not  used  for  these  equations.  Dr.  V.J.P.  Cullen  suggested  an 
ur.de-  relaxation  scheme  for  which  three  passes  over  the 
domain  should  produce  a  solution  of  acceptable  accuracy, 
since  the  coefficient  matrix  is  so  strongly  diagonally 
dominant.  .yass  lumping  of  the  coefficient  matrix  is  used  for 
the  first  guess.  This  technique  requires  replacing  the  mass 
mat-ix  [«]  by  the  identity  matrix  [I]  .  A  first  guess  of  this 
type  i?  *ble  to  describe  most  of  the  large  scale  features, 
whi~h  in  turn  reduces  the  number  of  iterative  passes  ever 
the  field.  Successive  passes  converge  to  solutions  which 
describe  smaller  scale  motion,  approximately  to  the  same 
order  of  magnitude  as  introduced  by  computational  error,  so 
that  further  iterations  are  not  needed. 
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I.   GRID  GIOPSTRT 

The  dor-air.  of  this  trodel  is  a  cylindrical  channel,  with 
total  length  of  4245  Km  and  width  of  3503  Km.  The  channel 
simulates  a  belt  around  the  earth  and  it  proves  to  be  an 
excellent  test  bed  for  comparing  with  the  finite  element 
formulations  used  by  Selley  (1576)  and  Older  (1981). 

The  domain  is  subdivided  into'  equilateral  triangles  as 
shown  in  Figure  6.  Post  of  the  test  runs  for  this  thesis  use 
a  12x12  mesh  which  has  156  nodes  and  289  elements.  This 
implementation  is  not  restricted  to  one  grid  pattern.  The 
technique  developed  by  Older  (1931)  to  vary  the  nodal 
geometry  smoothly  to  achieve  areas  of  denser  and  coarser 
resolution  is  also  implemented,  as  in  a  third  grid  pattern 
that  varies  the  nodal  geometry  abruptly.  A  short  discussion 
of  these  nodal  geometries  with  accompaning  illustrations  of 
each  is  presented  in  Chapter  VII,  where  the  different  test 
cases  are  described. 

Cyclic  continuity  is  assured  in  the  x  direction  by 
wrapping  the  domain  around  the  earth  to  form  a  cylindrical 
domain.  This  has  the  advantage  of  eliminating  the  east-west 
boundaries  and  it  simulates  the  flow  around  the  earth.  The 
only  boundaries  en  this  domain  are  the  north-south  walls  and 
their  treatment  will  be  discussed  shortly. 
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F.   INITIAL  CONDITIONS 

As  mentioned  previously,  the  reformulation  of  the 
governing    equations   into    the   vorticity-divergence 

shallow  water  equation  set  requires  solving  the  tirre 
dependent  variables  in  terms  of  the  stream  function  and 
velocity  potential.  The  continuity  equation  is  not  altered, 
so  that  its  solution  is  expressed  in  terms  of  i. 

For  the  basic  testing  of  the  model's  performance,  simple 
analytic  sinusoidal  initial  conditions  are  used  to  insure 
the  most  accurate  analysis  possible  and  to  simplify  the 
comparisons . 

The  sinusoidal  initial  fields  are  graphically  shown  in 
Figure  7  as  3-dimensional  surfaces.  The  geopotential  field  0 
consist  of  a  half  sine  wave  in  the  y  direction  and  a  single 
cosine  wave  in  the  x  direction.  The  stream  function  ¥, 
calculated  by  dividing  the  geopotential  field  by  the 
coriolis  force,  has  the  same  physical  structure  as  0.  The 
velocity  potential  1  has  a  single  sine  wave  in  the  x 
direction  and  a  half  sine  wave  in  the  y  direction. 

These  initial  conditions  are  computed  as  follows 

0  =   fQ  Asinoc1cosoc2  -  f0U(y  -  y^)  -  $ 

¥   =   0/fo  111-29 

1  =     CsinocjSinop  quasi-geostrophic   divergence 

where  A     -     arbitrary   amplitude 

&     -     coriolis    value   for   mid-channel    latitude 
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a)  0  and  Y  initial  fields. 


b)  X  initial  field. 
Figure        .  3~G-imer-sional  view  of  the  inital  fields. 
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3.   5CUNEARY  CONDITIONS 

3oundary  conditions  are  only  required  on  the  north  and 
south  walls  of  the  grid  domain.  Eue  to  cyclic  continuity, 
the  domain  is  wrapped  around  creating  a  cylinder  eliminating 
the  east  and  west  boundaries.  However,  careful  attention  to 
detail  is  needed  during  the  implementation  to  assure  this 
continuity.  Separate  boundary  conditions  are  applied  to  each 
of  the  predictor  equations  111-23,  111-24  and  111-25.  These 
conditions  are  computed  for  the  wall  nodes  only  and  are 
applied  during  each  pass  through  the  relaxation  scheme. 

The  vorticity  eauation  111-24,  the  most  sensitive  of  the 
predictor  equations  to  solve,  requires  ^  on  the  north-south 
boundaries  to  remain  constant  for  the  entire  forecast 
period.  Since  this  equation  is  solved  in  terms  of  ¥,  the 
Initial  north-south  ¥   values  are  saved  and  assigned   to   the 
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"boundary  points  after  each  pass  through  the  relaxation 
subroutine. 

-he  oroper  boundary  condition  for  the  divergence 
equation  III —25  would  be  dX/dn  =  C.  However,  for  the  purpose 
of  this  study,  there  is  more  interest  in  the  sinusoidal 
variation  in  the  y  direction  and  not  in  the  region  of  the 
walls.  Therefore  I  =  2   is  appropriate. 

The  continuity  equation  111-23,  the  most  complex 
predictor  equation,  requires  that  there  be  no  mass  flux 
through  the  north-south  walls.  The  geostrophic  boundary 
condition 

—   =   -uf0  111-30 

dy 

is   aoplied   to   the  north   south   boundary  nodes   for      the      terms 

[BEET]      in      equation   I 11-23 .    Integrating   the      inner      product 

<£  •■    V.,V. >   by   Darts   produces    the   boundary   terms 

If     v2(0.V    )V.    dxdy   =      f  /  *  •  (*  $  .V.  )  V  .dxdy 
yx  J    J  yx  J  J 

=    \[    [v.(Viv(^Vj))    -  v^V^-vV^    dxdy 

=      i  V.  v(0  .7.  )-S    ir   -    //   W.  .v(^.V,  )    dxdy 

=  r3ERTl  -  ^[<vjx.vix>  -  <Vjy.Vly>]  111-31 
where  n  is  a  unit  vector  normal  to  the  domain  and  dr  is  the 
differential  distance  along  the  path  of  integration  on  the 
perimeter   of    the  domain. 
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"he  geostrophic  boundary  condition  111-38  is  substituted 
into  the  contour  integral  in  equation  1 1 1-31  and  put  into 
laleriir  forrr,  in  the  sare  way  as  in  the  one-dimensional 
advective  equation  in  Chapter  II.  The  resulting  term  is 
derived  as  follows 


i      V.v'0.7.).n  dr   =  &   —  (#.V.)  V,  dx 


oy 


v*x 


(u  .  .  +  2u.  ♦  u.  .  ),   111-32 

3    J+l      J     J-l  i 

Equation  1 1 1-32  appears  twice  in  the  continuity  equation 
III  23,  for  *ime  levels  (n+1)  and  (n-1).  All  values  of  u  are 
known  "or  time  (n-1),  since  they  are  saved  from  the  previous 
calculations.  Fowever,  u(n+l)  has  not  been  computed.  To 
solve  for  u(n+l),  both¥(n+l)  and  T(n+1)  are  needed,  ^(n+l) 
is  solved  first  from  the  vorticity  equation.  1(n+l)  needs 
0(?»+l)  as  part  of  its  solution  and  0(n+l)  needs  u(n+l)  in 
its  solution.  To  avoid  this  problem,  it  is  assumed  that 
X(b+1)  has  a  negligible  contribution  to  the  solution  of 
ufn+1)  and  only  ¥( n*l )  is  used. 
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iv .     gcrruTiR  ieplimintation 

The  formulation  and  general  theory  of  the  finite  element 
method  was  presented  in  the  previous  chapters.  The  objective 
In  this  chapter  is  to  discuss  sere  Important  computational 
aspects  pertaining  to  the  implementation  of  the  finite 
element  prediction  system. 

The  nain  advantage  that  the  finite  element  method  has 
over  other  prediction  techniques  is  its  generality. 
Conceptually,  it  seems  possible  by  using  many  elements,  to 
approximate  virtually  any  surface  with  complex  boundaries 
iz1!  initial  conditions  to  such  a  degree  that  an  accurate 
solution  can  be  obtained.  In  practice,  however,  obvious 
engineering  limitations  arise,  a  most  important  one  being 
the  cost  of  the  computation.  As  the  number  of  elements  is 
increased,  a  larger  amount  of  computer  time  is  required  for 
a  forecast.  Furthermore,  the  limitations  of  the  program  and 
the  computer  may  prevent  the  use  of  a  large  number  of 
elements.  These  limitations  may  be  due  to  the  computer  speed 
and  stcrage  availability,  or  round-off  errors  propagated  in 
the  computations  because  of  finite  precision  arithmetic. 
Also,  the  malfunction  of  a  hardware  component,  if  the 
prediction  is  carried  out  using  many  computer  hours  to 
execute,  can  be  a  serious  problem.  It  is  therefore  desirable 
to  use  efficient  finite  element  programs. 
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The  effectiveness  of  a  program  depends  essentially  or. 
the  following  factors.  Firstly,  the  use  of  efficient  finite 
element  techniques  is  important.  Secondly,  efficient 
programming  methods  and  sophisticated  use  of  the  available 
computer  hardware  and  software  are  important.  The  third  very 
iirportfint  aspect  in  the  development  of  a  finite  element 
oro^ram  is  the  use  of  appropriate  numerical  techniques. 

The  vcrtici ty-divergence  model  described  in  the  previous 
chapter  is  implemented  on  the  13m"  3333  computer  located  at 
the  Naval  Postgraduate  School.  Some  notable  features  of  its 
architecture  are  the  three  trillion  bytes  of  virtual  mass 
storage,  of  which  eight  mega  bytes  are  available  to  each 
user,  and  the  c7  nanosecond  machine  cycle  time.  The  model  is 
executed  rostly  using  a  12x12  element  domain  requiring  4201c 
bytes  of  storage  and  30  seconds  of  CPU  time  to  execute. 
Irceeding  execution  time  and/or  available  storage  is  not  a 
problem,  in  fact  the  system  allowed  a  lot  of  flexibility 
during  the  implementation  phase  of  the  model. 

The  source  code  is  written  using  FORTRAN  IV  and  compiled 
on  an  optimizing  Fortran  H  compiler.  Appendix  A  contains  the 
source  cede  listing,  which  is  divided  into  five  subdivisions 
delineating  the  logical  structure  of  the  program. 

A.   PROGRAM  ARCEIT£CTURZ 

Program  features  incorporated  in  the  model  are: 

1)   Modularity.   With   only   a   few  exceptions,  each 
module  is  limited  to  one  page  of  FORTRAN  code.  This  makes  it 
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easier  to  comprehend  the  program.  Zach  module  performs  only 
o-^e  task.  For  example,  subroutine  CONTIC  computes  the  value 
jf  the  forcing  terms  for  the  continuity  equation.  likewise 
there  is  also  one  module  for  the  divergence  and  vcrticity 
equations.  To  implement  a  new  set  of  equations,  only  these 
modules  would  have  to  be  altered. 

"asily  controllable  switches.  Switches  may  be  set 
to  either  print,  plot  or  tabulate  harmonic  analysis  data  for 
most  of  the  available  fields.  The  ability  to  display 
intermediate  results  allows  each  portion  of  the  algorithm  to 
be  monitored  for  computational  adjustments.  This  also  makes 
it  easier  for  unfamiliar  users  to  become  acquinated  with  the 
computational  model. 

3 ^  Forcing  term  subroutines.  In  previous 
irpiementaticns •  each  forcing  term  was  calculated  by  a 
special  subroutine.  In  this  implementation,  the  calculations 
are  accomplished  by  general  purpose  routines  which  simplify 
the  implementation  of  the  complex  prognostic  equations:.  Hits; 
allows  implementation  of  different  equation  sets  (i.e. 
^aroclinic  yodel)  over  the  same  domain  with  minimal  effort. 

4<  T cementation .  Zach  variable  is  defined  by  a 
short  phrase  'Appendix  A,  A.).  The  function  of  each  module 
is  described  in  an  introductory  paragraph.  Shared  data  is 
placed  in  named  common  blocks  and  identified  with  each 
subroutine  which  uses  them.  A  subroutine  index  is  given. 


1 .   "air,  Program 

The  main  program  is  short,  calling  only  six  modules 
which  reflect  the  basic  sequential  flow  of  the  model.  It 
starts  with  initialization  of  all  model  parameters  (i.e. 
model  options,  domain,  finite  element  arrays,  inner 
products).  It  then  initializes  the  input  fields  (i.e. 
secpctential  heights,  stream  function  and  velocity 
potential)  and  is  followed  by  initialization  of  all 
remaining  dependent  variables.  At  this  time  the  model  is 
totally  initialized  and  time  integration  begins.  As 
mentioned  previously  two  techniques  are  employed  for  time 
integration,  each  having  its  own  module.  Upon  completion  of 
the  last  forecast,  the  program  terminates. 

Arrays  are  the  only  data  structures  used  and  are 
grouped  using  19  different  common  blocks.  Several  arrays  are 
use!  as  static  link  lists,  as  described  in  detail  later, 
which  simplified  the  algorithms.  The  common  block  format  has 
the  advantage  of  reducing  the  overall  execution  time  of  the 
program.  Most  of  the  arguments  passed  during  a  call  to  a 
subroutine  are  contained  in  conmon.  This  requires  less  time 
to  execute  since  no  parameter  passing  is  required  for  the 
arguments.  Another  benefit  of  this  format  is  that  the  code 
becomes  less  cumbersome  and  more  readable.  Each  variable  and 
array  is  defined  in  the  first  subsection  of  Appendix  A  along 
with  a  oage  index  for  the  subroutines. 


2.      Initialisation  Phase 

Appendix  A,  Section  C  contains  all  the  subroutines 
used  during  the  initialization  phase  of  the  program.  Frorr 
the  user  point  of  view,  the  most  important  suoroutine  is 
INI-S3,  the  first  subroutine  called,  which  contains  all  the 
global  variables  that  control  the  different  options 
available  per  run.  This  is  the  only  subroutine  that  is 
changed  to  run  the  different  experiments,  assuming  that  no 
new  computational  technique  is  introduced.  The  selection  of 
options  ere: 

1)  -  channel   location   -   the   channel  may  be 

shifted  north  or  south  by  presetting  the 
north/south  latitude  limits  in  INITG3. 

2)  -  variable  geometry  -  the  node  positions  may 

be  grouped  for  more  dense  node  patterns  to 
yield  higher  resolution.  Two  variables  HI 
and  H2  set  the  ratio  used  to  vary  the 
spacing  -along  the  x  .and  y.  axes, 
respectively. 

3)  -  initial  field  wave  length  and  amplitude  can 

be  altered  to  produce  various  effects. 

4)  change  the  initial  mean  flow. 

5)  -  diffusion   can    be   entered  for  any  of  the 

three  prognostic  equations. 

6)  -  maximum  length  of  forecast   period  may   be 

changed   and   a   print,   plot   or   harmonic 
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analysis  of  any  dependent  variable  may   be 
requested  for  any  tire  interval. 

Or.ce  the  experiment  is  determined,  the  options 
listed  above  are  set.  The  program  is  ready  to  be  executed. 

The  largest  part  of  the  initialization  phase 
consists  of  establishing  the  domain  and  producing  ail  the 
finite  elerent  computational  vectors  that  remain  constant 
throughout  the  experiment. 

"he  first  several  steps  in  setting  up  the  domain  are 
concerned  with  indexing.  Subroutine  COF.RES  is  called  first, 
where  all  the  nodes  (grid  points)  and  elements  (triangular 
areas'1  are  numbered  consecutively  starting  at  the  southwest 
corner  of  the  domain  and  moving  eastward  across  each  row  or 
latitude.  lor  earch  element,  a  record-  -of  all  of  its  nodes 
(▼ertlces)  *re  stored  in  array  ELMENT  (M,2),  where  M  is  the 
total  number  of  elements.  To  facilitate  the  inner  product 
evaluation  later,  a  local  numbering  system  is  required  for 
each  element.  That  is,  for  each  element,  its  nodes  are 
stored  counterclockwise  in  a  positive  sense.  The  first  node 
however,  is  arbitrary. 

With  the  domain  divided  and  numbered,  a  connectivity 
list  ( the  correspondence  between  each  node  and  the  neighbor 
nodes)  is  constructed  for  each  node  by  subroutine  CORRIL. 
lach  node  is  adjacent  to  four  or  six  other  nodes  depending 
on  whether  it  is  a  boundary  or  interior  node,  respectively. 
These   adjacent  nodes,  plus  itself,  make  up  the  connectivity 
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list  for  one  node.  The  connectivity  lists  are  then 
concatenated  sequentially  starting  with  the  first  nodes 
connectivity  list  into  the  vector  NAME  (NN),  where  NN  is  the 
sur  of  the  nodes  in  each  connectivity  list.  (i.e.  for  a 
12x12  dorain  with  156  nodes,  and  equilateral  elements  NN  = 
1044).  3or  the  first  tiire  during  the  initialization  phase, 
special  attention  is  given  to  cyclic  continuity.  As 
discussed  earlier,  cyclic  continuity  is  the  joining  of  the 
east  and  west  "boundaries  to  create  a  cylindrical  channel. 
The  connectivity  list  for  these  east/west  boundary  nodes 
must  be  complete  to  insure  proper  continuity  for  the 
calculations  later. 

The  connectivity  vector  NAM  is  frequently  used 
during  most  computations.  Two  utility  vectors  ISTART 
(containing  the  starting  location  in  NAM  for  a  particular 
node'  and  MUM  (containing  the  number  of  nodes  in  its 
connectivity  list)  are  used  to  locate  and  indei  through  the 
vector  NAMI,  as  will  be  seen  shortly.  This  same  technique  is 
used  to  index  through  the  coefficient  matrices  and  used 
during  most  of  the  node  interaction  computations. 

The  physical  properties  of  the  channel  are 
calculated  next  in  subroutine  CHANAL.  Here  the  north  and 
south  latitude  boundaries,  which  were  pre-set  in  INITG£  by 
the  user,  are  used  to  compute  the  grid  spacing  along  the  x 
and  y  axis.  Since  this  channel  simulates  a  belt  around  the 
earth,  the  magnitudes  of  both  LIITAX  and  DELTAY  (meters)  are 
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proportional  to  the  width  of  the  channel  divided  by  the 
number  of  ,?rid  points  in  the  y  axis. 

The  Cartesian  coordinates  for  each  node  are  computed 
by  subroutine  LCCATI  using  the  DILTAX  and  EILTAY  calculated 
in  CEANAL.  If  the  option  to  use  varying  grid  geometry  is 
i-sired,  subroutine  TRANS  transforms  the  grid  geometry. 
TRANS  also  computes  the  corresponding  new  Cartesian 
coordinate  values  for  each  prid  point  and  calculates  the 
minimum  E5ITAX  and  IELTAY  within  the  domain.  When  the 
secmet-y  is  changed  to  create  a  smaller  DELTAX  or  DELTA!, 
the  two  dimensional  advective  stability  criteria  is  also 
charged.  A  new  time  step  DT  has  to  be  computed  using 
equation  II 1-26.  Since  TRANS  transforms  the  geometry,  it 
also  computes  the  new  CT . 

Another  transformation  is  required  as  discussed  in 
Chapter  II.  The  transformation  from  Cartesian  coordinates  to 
area  coordinates  is  needed  to  perform  the  area  integration 
of  the  inner  products.  Sub-routine  AREA  computes  these 
transformations  exactly  as  outlined  in  Chapter  II,  Section 
C.  Again  cyclic  continuity  is  very  important  and  special 
care  is  needed  to  insure  proper  transformation. 

Following  the  area  transformation  is  the  computation 
of  all  the  inner  products  that  are  required  to  solve  the 
equations.  The  advantage  of  using  area  coordinates  is  that 
the  inner  products  (function  of  space  coordinates  only)  are 
computed   and   stored   once   and   used   repeatedly   without 
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recalculation-   Subroutine   INN2P.  computes  and  stores  these 
products  using  the  formulas  derived  In  Chapter  II. 

The  3oefficlent  matrix,  dimensioned  NxN  ,  where  N  is 
the  total  number  of  nodes,  is  a  matrix  of  coefficients  whose 
rows  are  the  equations  of  the  system  to  he  solved.  As 
discussed  In  the  computational  technique  section,  the 
members  of  this  sparse  matrix  are  the  inner  products  for  the 
left  hand  sides  of  the  equations.  Three  coefficient  matrices 
are  used  in  the  solution  of  the  equations.  The  diagnostic 
equations  (III-?  through  111-15)  use  a  coefficient  matrix 
with  the  inner  product  <V  ,,L  >  which  is  constructed  by 
subroutine  AMTHXl  and  stored  in  compacted  form  in  vector 
GfNN)  by  subroutine  ASEMBL.  Eowever,  when  solving  the 
prognostic  equations,  these  coefficient  matrices  have  a  DT 
ftlire  step  in  seconds)  term,  so  that  these  matrices  are  not 
assembled  until  the  time  integration  begins.  The  vorticity 
and  divergence  equations  (111-24,111-25)  use  the  coefficient 
matrix  H(N'N)  with  inner  products  <vjx.?lx>  +  <vjy*\y>  in 
solving  the  ?oisson  equations  for  the  stream  function  and 
velocity  potential,  respectively.  The  continuity  equation 
[111-23)  uses  a  combination  of  inner  products  in  its 
coefficient  matrix  J(N'N)  as  follows 


jx» v ix^   V¥jy  *  viy' 


J(At) 


<vj,vi> 


iv-i 
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to  solve  the  Helirholtz  equation.  At  the  start  of  each  tirre 
Integration  module,  subroutine  AMTHX2  is  called  to  construct 
the  two  coefficient  matrices  E  and  F. 

These  banded  and  sparse  rratrices  are  ccrrpacted  into 
vectors  to  save  storage  during  their  assemblage  by  ASEM3L. 
The  vectors  are  dimensioned  NN,  as  is  NAME,  the  connectivity 
vector,  and  both  use  ISTART  and  NUM  to  index  through  them. 
This  compaction  routine  was  used  by  Kelley  (1976)  and  Older 
'1981)  in  their  models,  but  was  developed  by  Hinsman  (1975). 

To  illustrate  matrix  assemblage  using  an  element   by 

element    technique,   consider   Figure  8.   Note   that   this 

illustration  is  for  element  number  3,  but  all   elements  are 

treated   in  a   similar  maner.   The  computational  technique 

required  that  for  each  point  (node)   describing  element   3, 

namely   nodes   Z,  3  and  14  stored  in  array  ELMENT,  the  inner 

nro-uct  <TT.,7J>  between  those  points  be  distributed  to  their 
J  i 

proper  location  in  the  coefficient  matrix. 

Subroutine  AMTRI1  builds  the  inner,  product  nodal 
interaction  and  stores  it  in  matrix  3,  dimensioned  3x3. 
Figure  6  illustrates  the  3  matrix  for  element  3,  where  the 
inner  product  <V  .  ,V.>  is  the  multiplicand  of  the 
corresponding  basis  and  test  functions,  respectively. 

The  local  dispensing  of  interactions  is  done  in 
ASZVBL.  Consider  the  second  row  of  [3]  in  Figure  8.  These 
are  the  interactions  between  node  3  of  element  3  to  the  test 
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AMTRX  -  builds  J3J  ^OT   one  element,  passes  [3]  and  the  element  number 
associated  with  [2]  to  ASEK3L.  This  process  continues  till  all  element 
node  interactions  are  assembled  in  the  coefficient  matrix. 
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ASEMBL  -  assembles  the  node  interactions  into  the  coefficient  matrix 

[Cj  ,  which  has  the  same  structure  as  HAKE.  The  following  diagram  assembles 

inner  product  '•▼jj. 


into  the  coefficient  matrix  [c] ,  for  element  J. 
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pseudo-code 

-do   i  -  1  -•  3     2 

11  -  SLMENT(l)  1  3 

u 

-do   j  -  START(ii)  -*  LAST(li) 

jj  -  NAME(j) 

-do   k  •  l-»3 

2 

kk  -  2LMEKT(k)  1 

3 

if  (  kk  -  jj  ) 

14 
then 

C(J)  -  G(j)  *   B<i,k) 

else  continue 

Lend  do 

-end  do 

-end 

do 
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Figure  8.  Assembling  and  storing  the  coefficient  matrix  for  element  3. 
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functions.  AS1KBI  locates  nodes  3's  connectivity  list  in 
NAM!   using   ISTARI   and   NUM.   In   Figure   8,   this  list  is 

delineated  ty  START(3)  and  LAST(3).  Now  ASEMBL  steps  through 
the  connectivity  list  for  three  iterations.  Turing  each 
pass.  ASI.V3L  is  searching  for  one  of  the  three  node  numbers 
for  element  3.  When  a  match  is  found  with  one  of  element  3's 
nodes  'i.e.  2,3  or  14)  and  node  3's  connectivity  list  (i.e. 
3,2,14,15  or  4)  the  proper  position,  to  which  this 
interaction  is  to  he  added  has  been  found  in  the  coefficient 
matrix.  Since  MAPE  and  vector  C,  the  compacted  coefficient 
rat  Mi,  are  dimensioned  identically,  the  same  pointer  (i.e. 
J  in  Figure  3)  is  used  to  index  through  both  arrays.  This 
procedure  is  repeated  for  all  elements  in  the  domain  to 
assemble  the  coefficient  matrix  of  the  equations.  The 
pseudo-code  for  ASEMBL  is  shewn  in  Figure  8  to  facilitate 
stepping  through  this  exarrple. 

The  domain  and  all  finite  element  work:  vectors  are 
initialized  at  this  point.  Subroutine  ER-MS3T  is  called- .later 
to  compute  interpolation  points  for  the  harmonic  analysis 
subroutines. 

The  last  phase  of  the  initialization  process  is  the 
initialization  of  the  dependent  variables.  The  three  input 
fields  geopotentlal  heights,  stream  function  and  velocity 
potential  are  computed  in  subroutine  IC  using  the  equation 
set  II 1-30 .  However,  the  variables  calculated  from  the 
diagnostic   equations   have   to   be  computed  using  the  input 


fields.  These  varieties  are  used  during  each  tiire  step  while 
solving  the  prognostic  equations. 

The  diagnostic  equations  are  solved  in  subroutine 
DIPTAB,  first  during  the  initialization  phase  and  later 
during  the  time  integration  phase.  Each  diagnostic  equation 
calls  its  own  module  to  compute  the  value  of  the  forcing 
function  and  stores  the  computed  values  in  the  vector  RES. 
These  equations  all  use  the  same  coefficient  matrix  when 
solving  the  diagnostic  equations.  Subroutine  SOLVER  is 
sufficiently  genereal  to  solve  each  equation.  SOLVER  uses 
vector  RES  and  coefficient  matrix  G  to  under-relax  towards 
the  solution.  As  mentioned  previously,  the  coefficient 
matrix  is  strongly  diagonally  dominant  so  that  three  passes 
over  the  domain  are  sufficient.  At  the  end  of  rEFVAR,  output 
is  generated  depending  on  what  print,  plct,  or  harmonic 
analysis  controls  were  preset. 

This  completes  the  initialization  pnase  of  the  model 
and  the  program,  for  the  forecast  phase  will  he  described 
next . 

2 .   Forecast  Phase 

The  forecast  phase  is  accomplished  in  two  steps.  The 
first  time  step  is  made  using  two  half  steps  by  subroutine 
r*£TZNO.  Eere  the  prognostic  coefficient  matrices  are 
constructed  using  half  the  ET  value  by  calling  AMTRX2. 
Ar"TFX~  uses  the  same  computational  technique  to  construct 
the  ccefficient  matrices  as  described  for  AMTRX1 . 
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lach  of  the  prognostic  eauations  111-23,  111-24  and 
111-25  calls  its  own  subroutine  (CONTIC,  VORTEQ  and  TIVEC 
respectively'  to  compute  all  the  terms  on  the  right  hand 
side,  which  are  stored  in  the  vector  HHS.  After  computations 
for  RES  are  completed,  subroutine  RELAX  solves  the  equations 
by  over-relaxation  as  described  in  the  computational 
technique  in  Chapter  3.  Once  the  solutions  for  the  (n+1) 
time  step  are  completed,  DEPAB  is  called  to  update  the 
variables  from  the  diagnostic  equations.  Two  oasses  through 
rtATZNO  advances  the  solution  fields  one  time  step. 

The  remainder  of  the  forecast  period  is  integrated 
using  the  leapfrog  scheme.  Subroutine  L2APFR  performs  this 
integration  using  the  identical  format  as  vATZNO,  except 
that  IT  equals  two  DT.  At  preset  times  as  specified  in 
INITG5,  the  different  fields  are  saved  for  printing.  This 
process  continues  until  the  final  forecast  time  is  reached. 

B.   UTILITY  YOTULZS 

Ir.ce  the  equation  formulation  is  completed,  as  in 
Chapter  III,  all  the  inner  products  and  types  of 
integrations  are  inown.  Versatile  modules  can  oe  written  to 
per^crm  these  computations.  Consider  a  term  of  general  form 
<^A  .  V  .  ,V  .  >  where  i  is  the  node  about  which  the  term  is 
evaluated  and  the  j's  are  the  nodes  connected  to  node  i,  or 
the  surrounding  nodes.  The  inner  product  values  <V.,V  >  are 
already  computed  and  stored  for  all  the  nodes,  during  the 
Initialization  nhase  of  the  model. 
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-he  regaining  computation  to  complete  trie  evaluation  of 
this  term  is  the  mult  iplicat  ion  of  the  scalar  coefficient  of 
A  it  node  1  with  the  corresponding  inner  product  <V.,V. >  for 
node  «.  This  requires  indexing  through  node  i 's  connectivity 
list  stored  in  vector  NAME,  and  for  each  node  in  the  list 
irnltiply  and  total  the  products.  The  cumulative  sum  of  these 
multiplications  is  assigned  as  node  i 's  contribution  for 
this  term.  Subroutine  TERV3  performs  this  exact  computation. 
All  that  is  passed  to  7ERM3  is  the  scalar  field  A  and  the 
sign  of  the  inner  product,  TEBM3  then  computes  the 
contribution  for  each  node  in  the  domain  and  accumulates  it 
in  the  wcric  vector  RES. 

Three   other  utility  modules  are;  TERM1,  which  computes 

the  first  scalar  multiplication  for   triple   inner   products 

'i.e.  /Aj  Vj3kVk  ,Vi>) .  The  product  <V.Vk,V£f  is  again  already 

computed  and  stored  by  subroutine  INNER.  TERM1  computes  <\ ?. 

Vk  ,7 i >   and   constructs   a  compacted  vector  similar  to  the 

coefficient  matrices.  This  reduces  the  effort  of  multiplying 

the  second  scalar  to  a   TERf3   computation.   TZRP2   computes 

node  interaction  of  the  following  type  <Aj  Vjxt 7^x> ,  where  both 

the   basis   and   the  test  functions  are  derivatives.  Lastly, 

subroutine  TIRK4  computes  node  interaction  for  terms  as  <A.V. 

3    Jx 

,v\  ,  where  only  the  basis  function  is  a  derivative. 

When  examining  the  right  hand  side  of  the  equation  sets 
111-23,  111-24  and  111-25,  it  is  obvious  this  implementation 
is  a  subscripting  nightmare;  however,  the  use  of  the  utility 
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modules  TESM1,  TERM2 ,  TERM3  and  TERM4  simplifies  the 
irplementation  to  only  determining  what  order  to  call  the 
utility  modules.  Examination  of  subroutines  CONTEQ,  VCRTEQ 
and  EIVEC,  which  compute  the  right  hand  sides  for  111-23, 
111-24  and  1 1 1-25  respectively,  illustrates  this  fact.  No 
other  subroutines  or  calculations  were  required. 
Implementation  of  these  equations  required  minimal  effort. 
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V.   P3IMITIYI  MOBIL  EXPERIMENT 

The  previous  two  chapters  presented  the  detailed 
formulation  and  implementation  of  the  vorticity-di vergence , 
shallow  water  equation  rrodel.  The  results  frorr  this  rrodel 
will  be  compared  with  the  results  from  the  primitive  model 
in  Charter  711.  To  facilitate  interpretation  of  the 
comparisons,  a  "brief  description  of  the  primitive  model 
follows.  See  lelley  (1976)  for  a  detailed  discussion  of  the 
entire  model . 

Also  presented  in  this  chapter  is  an  experiment  which 
demonstrates  significant  improvement  of  the  solution  from 
the  primitive  model.  lelley's  implementation  used  elements 
which  were  right  triangles.  Cider  (1981)  showed  that 
eauilateral  elements  are  far  superior  to  triangular 
elements.  This  experiment  re  implements  the  primitive  model 
using  equilateral  elements  and  a  comparison  is  made  between 
the  results  cf  both  implementations. 

A.   ,VC-IL  DESCRIPTION 

A  form  of  the  barctropic,  shallow  water,  primitive 
equations  developed  by  Phillips  (1959)  is  used  as  the 
£Overnin£  equation  set  for  this  model.  In  Cartesian 
coordinates  the  equation  set  is 


68 


—  = (u#) (v0) 

3t      dx         dy 


V-l 


3u  d0  du     du 

—  =  -  —  -u  —  -v  —  «■  vf0 

3t  ox  dx      c-y 

3v  d2f  dv     dv 

^t  dy  dx      dy 


V-2 


V-3 


The  finite  element  formulation  of  this  set  of  equations 
evaluates  the  height  and  velocity  components  at  the  same 
nodal  points.  This  is  an  important  consideration,  "because 
the  other  models  ( the  linearized  model,  see  Chapter  VI,  and 
the  vorticity-divergence  model,  see  Chapter  III)  either 
stagger  the  dependent  variables  or  have  the  property  of  a 
staggered  formulation.  When  comparisons  are  made  between  the 
models,  it  is  this  lattice  structure  that  is  being  compared. 

This  form  of  the  shallow  water  equations  includes 
gravity  waves  as  a  solution.  Gravity  waves  have  a  maximum 
phase  soeed  of  about  300  meters/second.  When  the  correct 
tine  step  is  calculated  using  equation  111-26,  a 
considerably  smaller  time  step  is  obtained  compared  to  the 
larger  time  step  permitted  in  the  vorticity-divergence 
formulation.  This  is  an  important  feature.  If  solutions  from 
all  models  are  equally  as  good,  the  best  formulation  would 
be  determined  using  the  computational  time  required  to 
produce  the  desired  forecast. 


69 


All  Todels  use  the  same  domain  structure.  In  fact,  the 
domain  described  in  Chapter  III  was  patterned  after  the 
iomain  implemented  by  Kelley.  Again,  this  domain  simulates  a 
belt  around  the  earth,  with  cyclic  continuity  which 
eliminates  the  east-west  boundaries.  Rigid  boundary  vails 
exist  at  the  equator  and  at  30  degrees  north  latitude.  The 
domain  is  composed  of  a  12x12  point  mesh  and  subdivided  into 
the  right  triangular  elements  illustrated  in  Figure  9. 
Notice  that  the  grid  points  are  not  shifted  as  in  the 
equilateral  element  implementation  shown  in  Figure  6. 

The  following  boundary  conditions  are  imposed: 

1)  -  no    cross   channel   flow   at   the  latitude 

boundaries . 

2)  -  a   geostrophic   balance   at   the   channel  walls 

imposed  on  the  continuity  equation  V-l. 

This  model  has  a  simple  second  order  diffusive  term  in 
the  equations  of  motion  V-2  and  V-3.  However,  for  the 
purpose  of  evaluating  these  different  formulations,  this 
option  was  not  implemented  during  the  comparison  phase. 

Initial  conditions  consist  of  a  single  wave  in  the  x 
direction  and  a  half  wave  in  the  y  direction.  The  initial 
fields  for  the  three  dependent  variables  are  shown  in  Figure 
12.  The  maxixum  zonal  wind  perturbation  of  5.5  meters/second 
is  superimposed  on  a  mean  zonal  flow  of  12   meters/second. 
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Figure  10.  Initial  fields  for  the  primitive  model.  Both  the  x  and  y 
axes  are  multiplied  "oy  10  Km. 
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This  cursory  description  of  the  primitive  model  mentions 
only  those  significant  features  that  will  weigh  heavily  in 
the  ccirparisions  later.  The  Galerkin  implementation  of  this 
model  is  similar  to  that  presented  in  Chapters  II  and  III 
and  the  system  of  equations  is  solved  using  a  Gauss-Seidel 
iterative  procedure.  Further  details  concerning  this 
primitive  model  are  given  in  Kelley  (1976). 

2.   RZSU1TS 

This  experiment  involves  the  shape  of  the  elements.  As 
mentioned  previously,  Kelley's  implementation  subdivides  the 
domain  into  right  triangular  elements,  as  illustrated  in 
Figure  9.  Considerable  small-scale  noise  was  observed  by 
Kelley  in  fhe  43  hour  forecast  solution. 

The  transition  from  right  triangles  tc  equilateral 
triangles  changes  the  size  of  the  domain.  With  right 
triangles,  the  ax  arday  grid  spacings  are  equal  (300  KM).  A 
12x12  grid  matrix  has  a  length  and  width  of  3600  KM.  With 
equilateral  triangles,  the*x  and  ay  grid  spacings  are  no 
longer  equal.  Arbitrarily,  the  *y  grid  spacing  is  held 
constant  (300  KM)  and  a  new^x  grid  spacing  computed  by 

u  =   Ay/ cos (30)  V-4 

A  12x12  grid  matrix  with  equilateral  elements  has  a  width  of 
3600  KM  and  a  length  of  4045  KM. 

Figure  11  contains  the  48  hour  forecasts  produced  using 
both  types  of  elements.  The  three  dependent  variables  fields 
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Figure  11.  US   hour  forecast  comparison  from  the  primitive  model  using 
both  right  triangles  and  equilateral  triangles  for  element  shapes. 
Arbitrary  perturbation  amplitude  of  5*5  m/s. 
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are  compared  for  each.  The  small -scale  noise  that  Xelley 
observed  is  present.  The  three  fields  shew  varying  degrees 
cf  distortions,  vhich  are  especially  noticeable  along  the 
"boundaries  . 

Cider  (1981)  found  that  theroot  rrean  square  error  (RMSE) 
was  reduced  22  percent  by  using  equilateral  shaped  elements. 
This  improvement  is  apparent  on  viewing  Figures  lib,  d  and 
f.  The  contours  are  smooth  and  the  boundaries  are 
noise  free.  Eelley  showed  excellent  treatment  of  wave 
propagation  by  this  primitive  model.  The  lowest  resolution 
grid  '6x6)  tested  by  Kelley  was  within  four  percent  of  the 
true  phase  velocity.  Changing  the  element  shape  had  no 
apparent  effect  on  the  phase  velocities. 

recause  the  outcome  cf  this  experiment  was  a 
significantly  improved  forecast  solution,  future  comparisons 
with  the  primitive  model  will  be  made  using  equilateral 
elements . 
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"":  e  previous  chapter  demonst rated  how  the  shape  of  the 
element  can  significantly  improve  the  solution.  Williams  and 
?ienktewicz  (1951}  used  differently  shaped  basis  functions 
c  -  a  linearized  e o ua t i o n  set  to  o r o d uc e  excellent  solutions 
when  applied  to  the  gecstrophic  adjustment  problem. 

Tra'ial  staggering  of  dependent  variables  in  finite 
difference  formulations  has  given  rucn  better  solutions  to 
*'re  geostrophic  adjustment  orocess,  and  these  forms  are 
widely  used  in  meteorology.  Schoenstadt  (1930)  found  similar 
result?  with  finite  element  formulations  with  piece wise 
linear  basis  functions.  However,  staggering  nodal  points  is 
net  a  convenient  rethod  to  implement,  especially  in 
two-dimensions  with  irregular  boundaries*  so  alternative 
schemes  are  needed  . 

The  implementation  of  the  alternative  scheme  introduced 
by  'illia-s  and  Zienkiewicz  (1931!  are  presented  in  this 
chapter.  As  mentioned  above,  this  formulation  uses  different 
basis  functions  for  the  height  and  the  velocity  fields.  One 
of  rhe  basis  functions  is  piecewise  linear,  while  the   other 


is   piecewise  constant,  as  is  illustrated  in  ? i g u r 
ore  ^TO'ciT^1  do  ma  in  • 


S   1  ~ 


for  a 


Yfc 


1  - 


a)  Piecewise  linear  basis  function. 


-I h 


^-AX->I 


b)  Piecewise  constant  basis  function 


Figure  12  .   Different  shaped  basis  functions. 
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A.   ICUATION  REFORMULATION 

The  priritive  forrr  of  the  shallow  water  equations 
^resented  in  Chapter  V  (V-l.V-2  and  V-3)  is  used  to  derive 
the  linearized  equations  needed  for  this  experiment.  The 
velocity  equations  V-2  and  V-3  reirain  unchanged  and  a  linear 
basis  function  (V  .  )  is  used  to  approximate  the  u  and  v 
variables . 

The  continuity  equation  7-1  is  linearized  as  follows: 


VI-1 


where  f  is  the  average  geopotential  over  the  domain.  A 
plecewise   constant    basis    function   (V  .  )   is  used   to 

u 

approximate  the  geopotential.  This  linearization  is 
reasonable  in  this  case  because  the  Rossby  radius  of 
Reformation  jj  f\  is  much  larger  than  ax  [see  Williams  and 
Zie^kievlez  (1981)].  The  Galerkin  method:  is  applied  to  this 
linearized  equation  set  using  a  plecewise  linear  test 
function  for  V-2  and  V-3,  and  a  plecewise  constant  test 
function  for  VI-1. 

The  oiecewise  constant  basis  function  has  the  property 
cf  displacing  the  geopotential  to  the  centroid  location  of 
the  elements,  which  should  give  the  same  effect  as 
staggering  the  grid  points,  as  seen  in  Figure  13.  The 
density  of  geopotential  data  in  the  domain  is  now  greater 
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Figure  13  .  Staggering  of  the  geopotential 
about  the  velocity.  The  a   symbol  are  the 
actual  grid  point  locations  and  the  •  symbol 
the  centroid  location  of  each  element. 
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than   the  density  of   the   velocity  data.   Instead  of   a 
gecoctent  ial   value   for   each  node,   there  is  one  averaged 


val*.te  over  each  element. 


""he  final  form  of  the  Galerkln  equation   for   VI-1,   7-2 
aid  7-3  after  performing  the  tire  differencing  is: 


**l,  v      -  -    ..n-1 


u  ;  <v4f7 


i>  -  u     <vj-';i>  -  •*  »$<**, -V  ♦  *Kay"*x'h> 


+  v.u ,  vV  .  V.    V  ">  -  f  vn<V   V  >1 


VI  -3 


r^/^ 


n  n  .  „ 


V    V  > 
71-4 


This  linearized  set  of  equations  71-2,  71-3  and  71-4  is 
solved  using  a  ^auss  Seidel  iterative  procedure.  It  is  worth 
mentioning  that  the  coefficient  matrix  <W. ,V.>  in  equation 
7!-2  has  all  non  zero  coefficients  equal  to  one,  since  the 
integration  of  the  inner  product  <W.,Wi>  involves  piecevise 
constant  functions. 

This  equation  set  is  implemented  rather  than  the 
equation  set  7-1,  7-2  ar.i  7-3  using  all  of  the  existing 
primitive  model  cede.  The  major  modification  involved  the 
way  that  the  geopotential  was  referenced.  With  an  average 
geopotential  over  each  element  instead  of  a  value  at  each 
node  fcr  a  12x12  mesh,  there  are  233  geopotential  points 
versus  the  156  velocity  (u,v)  points. 


E.   RESULTS 

The  results  from  this  linearized  model  are  compared  to 
those  from  the  primitive  model.  The  initial  field  for  this 
experiment,  Figure  14  has  a  maximium  perturbation  zonal  wind 
that  is  one  fifth  of  the  value  used  in  Chapter  V,  Figure  10 . 
The  mean  zonal  wind  remains  10  meters/second.  The  48  hour 
forecast  solutions  for  each  field  are  compared  in  Figure  15. 

This  alternative  formulation  shows  some  promise, 
although  there  are  seme  minor  perturbations  in  these 
contours  compared  to  those  in  the  primitive  solution.  No 
explanation  is  offered  for  this  small-scale  noise,  although 
possibly  the  other  formulations  presented  b/  Williams  and 
Zlenkievlcz  (1981)  would  improve  the  solutions. 
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Figure  1A.   Initial  fields  for  both  the  primitive  model  and  the 
linearized  model. 
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.:s  in  the  previous  two  chapters,  a  comparison  between 
two  models  will  be  presented.  The  results  from  the 
vort  ici  ty-di  vergence  model  are  ccmnared  to  the  results  fror1 
the  primitive  model.  To  fully  exploit  the  differences 
between  the  models  performances  and  indicate  the  strengths 
and  weaknesses  of  each  formulation,  three  domain  geometries 
i?A  three  initial  conditions  are  used.  All  solutions  are  at 
4S  hour,  except  for  one  case  which  was  extended  to  96  hours. 

Zrom  these  two  finite  element  formulations  some 
additional  insight  is  obtained  concerning  the  execution  time 
required  as  the  grid  resolution  changes.  Lastly,  a  "brief 
discussion  on  the  sensitivity  of  the  computational  technique 
is  p i ▼ e n . 

A.   TZST  DOMAINS  ANl  INITIAL  CONEITIONS 

The  three  domain  geometries  used  in  the  model  evaluation 
are  illustrated  in  Figure  16.  All  domains  consist  of  a  12x12 
element  -resh  with  equilateral  shaped  elements  (156  grid 
points  and  cyclic  continuity  is  imposed  on  the  east  and 
vest  boundaries.  The  domain  has  dimensions  of  4045  KM  along 
the  x  axis  and  35^3  KM  along  the  y  axis. 

The  regular  domain  'Figure  16a)  has  a  uniform 
distribution  of  -rid  ooints,  with  a  minimum  grid   spacing 
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Figure   16 .      Test  domain  geometries . 
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along  the  x  axis  of  337  £M .  This  is  the  most  congenial 
loiraia  ar.d  produces  the  best  results. 

The  smooth  domain  [Figure  16b)  has  a  smoothly  varying 
iistrioution  of  grid  points.  This  technique  allows  a  smooth 
variation  of  resolution.  It  was  developed  by  Cider  (1981), 
who  showed  how  it  significantly  reduced  noise  at  all 
frequencies  oo^pared  with  other  variable  grid  domains. 

The  degree  of  resolution  variation  is  accomplished  by 
selecting  an  appropriate  value  for  the  ratio  of  maximum 
Stretch  to  minimum  shrink:  along  both  axes.  The  experiments 
presented  in  this  chapter  use  a  2.5  ratio  for  both 
directions.  This  produces  a  grid  point  concentration  in  the 
right  -enter  of  Figure  16b  and  coarse  resolution  at  the  top 
srd  bottom  left  of  the  domain.  The  minimum  grid  spacing 
along  the  x  axis  is  159  Km:  in  the  dense  grid  point  area. 

The  third  domain  was  the  least  hospitable  geometry  for 
both  models.  The  abrupt  domain  (Figure  16c)  has  a  dense  grid 
point  concentration  on  the  left  and  coarse  spacing  en  the 
right  of  the  domain.  The  minimum  grid  spacing  along  the  j 
axis  is  15S  KM  in  the  area  on  the  left.  The  grid  spacings 
are  uniform  except  for  the  abrupt  change  along  the  center. 

Although  these  three  domains  are  simple  in  structure, 
they  are  adequate  to  test  both  models'  performance.  To 
further  enhance  some  differentiating  characteristics  between 
the  two  formulations,  three  initial  conditions  are  used.  Two 
have  previously   been  described  in  Chapters  V  and  VI.  Their 
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rain  iistinguishing  feature  is  the  arbitrary  perturbation 
airolitu<?e  ■  AVA)  magnitude.  The  nearly  linear  case  Las  an  APA 
=  1.1  m/s  compared  tc  5.5  m/s  for  the  more  nonlinear  case. 
All  initial  conditions  have  a  12  sn/s  mean  flew  component. 

tfith  the  nearly  linear  case  both  models  behave  well.  The 
introduction  of  the  rrore  nonlinear  initial  disturbance 
illustrated  the  boundary  and  computational  technique 
sensitivity.  The  third  initial  condition  is  the  nearly 
linear  initial  field  with  the  wave  length  equal  to  half  the 
domain  length.  This  has  the  effect  of  producing  two  waves 
with  the  domain  . 

2.   ^IST  CASI  COMPARISONS 

The  comparisons  between  models  are  divided  according  to 
the  domain  geometry.  The  primitive  model  has  three  fields! 
geopotential ,  a  and  v.  The  vorticity  -  divergence  model  has 
seven  fields:  geopotential,  stream  function,  velocity 
potential,  u,  v,  vorticity  and  divergence.  The  geopotential, 
u  ^tA    v  will  be  the  only  fields  used  for  the  comparison. 

1 .   Regular  Case 

The  regular  domain  geometry  (see  Figure  16a)  is  used 
for  this  first  comparison.  The  initial  conditions  have  an 
A?A  =  1.1  m/s,  as  shown  in  the  contour  plots  in  Figure  17. 

The  4:9  hour  solutions  for  both  models  are  presented 
in   Figure   16.   As   anticipated,   all  of  the  solutions  have 
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smoothly  "ontoured  fields  with  no  noticeble  noise  and   their 
phase  velocities  are  comparable. 

v-th  the  nearly  identical  solutions,  the 
distinguishing  feature  between  the  models  is  the 
computational  tire.  The  computational  time  is  determined  by 
the  size  of  the  time  step  each  model  is  allowed  to  use  as 
indicated  on  equation  111-26.  For  this  domain  the  gri i 
snacin*  is  33?  £* .  The  maximum  phase  velocity  possible  is 
different  for  each  model.  The  primitive  model  allows  gravity 
waves  so  c  -  300  m/s  whereas  the  vorticity-di vergence  model 
filters  out  these  high  frequency  waves.  A  c  =  12  m/s  is  used 
for  this  formulation. 


^odel 
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?"] 


Table  1 

<3X 

c 
(m/s) 

At 
(sec) 

it  7 
-z-77 

3?0 
10 

455 
13,755 

*   of  steps 
(iterations) 

37c 

13 


CPU  units 
(  sec) 


145. 


-  p. 


Coinoarisons  of  the  computational  times  for  a  4B  hour 
forecast  between  the  primitive  model  and  the 
vortinity-divergence  model. 


Table  1  comoares  the  results  of  the  computational 
times  for  both  models.  The  vorticity  divergence  model 
produced  as  accurate  results  over  the  regular  domain  using 
22*-  of  the  computational  time  needed  by  the  primitive  model. 
In  fa^t  from  geostrcphic  reasoning  the  vort icity-divergence 
model  should  produce  better  forecasts  when  small  scale 
features  are  oresent. 


2.   Smooth  Case 


c  rp  -\  r\  t  V 


omain  geometry  (Figure  16c 


s  used  in 


this  second  comparison.  The  initial  fields  shown   in  Ji^rure 
IS  have  a  disturtar.ee  amplitude  of  1.1  tr/s. 

The  48  hour  solution  comparison  is  presented  in 
Ji^'jre  22.  All  fields  again  have  srrooth  contoured  plots. 
Close  inspection  of  the  two  u  fields,  Figures  20c  and  i. 
shows  that  the  urimitive  u  field  has  small  kinks  along  sore 
contours  and  a  weak  tilt  near  the  central  boundary  ncdes, 
although  this  may  he  a  function  of  the  plotting  routine. 
Notice  the  £ood  symmetry  for  the  vort ic ity-di vergence  u 
field. 

As  in  the  regular  case,  this  smooth  experiment 
produced  twe  acceptable  solutions.  Again,  the  computational 
time  is  the  differentiating  criteria. 


Table  2 

vodel 

AX 

[EM) 
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(T/S) 

At 

(sec) 
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12 

271 

£124: 

n   of  steps 
(iterations ) 

c38 
21 


CPU  uj 

nits 

(  sei 

c) 

324 

.3 

49 

Comparison  of  the  computational  times  for  a  43  hour  forecast 
"etwee",  the  primitive  model  and  the  vort  icity-divergence 
model  using  the  smooth  domain. 


Table  2   ccm-oares   the  computational  times  for  both  models. 
The  vcrtici tv-diver*ence  has  a  33.8^  saving  of  CPU  time. 
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Figure  19.   Initial  fields  for  both  the  primitive  and  vorticity- 
ilvergence  models  using  the  smooth  domain  and  perturbation  amplitude 
of  1.1  m/s. 
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In  an  effort  to  contrast  the  computational  accuracy 
at  the  rodels,  this  experiment  is  repeated  using  the  rore 
nonlinear  initial  case.  Figure  21  shews  the  initial  fields 
with  a?A  =  5.5  rr/s .  This  larger  initial  disturbance  is 
reelected  in  the  greater  geopotential  amplitude  and 
magnitudes  of  the  contour  lines. 

Figure  22  shows  the  45  hour  solution  comparison.  As 
in  the  rore  linear  case  presented  above,  all  of  the  plots 
have  smooth  contours  with  no  noticeable  noise.  The 
vorticity-di vergence  geopotential  field,  Figure  22b,  has  the 
rid^e  extending  farther  north,  and  flattening  of  the 
southernmost  contour,  than  does  the  primitive  geopotential 
field,  Figure  22a.  The  mean  geopotential  heights  for  both 
models  have  also  increased.  The  primitive  mean  geopotential 
is  now  49£8<?  gum  and  the  vortici  tv-divergence  geopotential 
is  49630  gpm. 

These  two  discrepancies  indicate  that  the  boundaries 
are  not  handled  accurately.  Turing  the  implementation  of  the 
vortici ty-di vergence  model,  treatment  of  the  boundaries  was 
the  most  troublesome  phase.  The  vortici ty-divergence 
formulation  is  a  complex  equation  set  and  time  limitations 
restricted  further  investigation  of  mere  sophisticated 
boundary  conditions. 

This  same  initial  condition  is  now  extended  to  a  96 
hour  forecast,  which  is  shown  in  Figure  23.  The  mean 
vcrticity-divergence   2-eopotentiai   increased   to  49S00  gprr , 
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Figure  21 .  Initial  fields  for  both  the  primitive  and  the  vorticity- 
divergence  models  using  the  smooth  domain  and  perturbation  amplitude 
of  5.5  m/s. 
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Figure  23  •  96  hour  forecast  comparison  "between  the  primitive  model 
and  the  vorticity-divergence  model  using  the  smooth  domain.  (APA  = 
5-5   m/s) 
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whereas  the  mean  primitive  geooo tential   regained   constant. 
All   fields  have  smooth  contours.  Close  inspection  of  both  u 


?  * 


e  ..  z 


A  —  -4  .   W 


'2c   and  * 


i,  shows  a  slight   skewing  of   the 


contours  along  the  central  channel  grid  points.  The 
vorticity  divergence  u  field   has   the   rost   pronounced 

deviations.  A  hypothesis  for  this  skewing,  explained  further 
below,  is  that  it  is  caused  by  the  computational  technique 
employed.  The  relaxation  scherres  are  extremely  sensitive  and 
fine  tuning  of  the  relaxation  coefficient  would  nave 
rehired  more  time  than  was  available. 


'able  3 


Vodel 

A  x 
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(m/s) 

4  t 

(sec) 

199 
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2Z2 

12 

271 
8124 

#  of  steps 
( iterations ) 

1276 

42 


CPU  units 
(sec) 

594.2 
91.0 


Comparison  of  the  computational  times  for  a  9€  hour  forecast 
between  the  primitive  model  and  the  vortici ty-divergence 
Tccel  using  the  smooth  domain. 


Table  3  shows  the  comparison  between  troth  models' for 
the  96  hour  forecast.  A  savings  of  35  percent  is  realized 
with  the  vorticity  divergence  model. 

The  above  experiment  points  out  the  two  areas  where 
the  vorticity — divergence  model  is  presently  weak,  the 
Increase  of  the  geopotential  and  the  sensitivity  :f  the 
relaxation  coefficients.  3oth  these  weaknesses  can  be 
improved  and  are  not  a  result  of  the  formulation  but  of  the 
implementation .  At  present,  their  influence  is  not   detected 
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ex-ep'  In  extremely  long  forecasts,  such  as  in  the  96  hour 
exaxpie.  Further  experiments  ray  still  be  required  if  they 
demonstrate  significant  differences  in  the  solutions  cf  the 
two  models. 

The  last  experiment  on  the  smooth  domain  uses  the 
mere  linear  initial  condition  AFA  =1.1  m/s,  tut  its  wave 
length  is  divided  by  two,  so  that  two  shorter  waves  are 
propagated  through  the  channel.  The  initial  conditions  are 
shown  in  Figure  24.  Decreasing  the  wave  length  has  the  same 
effect  as  decreasing  the  density  of  grid  points.  In  this 
case  six  grid  points  are  used  to  describe  the  wave  structure 
versus  the  12  used  in  the  previous  cases. 

The  46  hour  comparisons  are  shown  in  Figure  25.  As 
in  all  previous  cases,  a  computational  time  saving  of  S4%  is 
rained  with  the  vcrticity-divergence  model.  With  fewer  grid 
points  describing  the  wave  structure,  more  small-scale  noise 
is  introduced  into  the  solution.  Comparing  the  primitive 
geopotential  field.  Figure  25a,  tc  the  initial  geopotential , 
Figure  24a,  shows  a  dampening  cf  the  wave  amplitude,  whereas 
the  vo rtici ty-divergence  geopotential.  Figure  25b, 
correlated  well  with  the  initial  geopotential. 

The  high  frequency  noise  is  evident  on  both  u 
fields,  Figures  25c  and  d.  The  primitive  model  u  field  is 
poorly  defined  along  the  boundary  and  becomes  irregular  over 
the  interior  grid  points. 
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Figure  2fr.   Initial  fields  for  both  the  primitive  and  the  vorticity- 
divergence  models  using  the  smooth  domain.  There  are  two  waves  embeded 
in  the  flow  and  A?  A  =  1.1  m/s. 
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Figure  25.  ^  hour  forecast  comparison  "between  the  primitive  model 
and  the  vorticity-divergence  model  using  the  smooth  domain.  Two  waves 
are  em ceded  in  the  flow  and  A?A  =  1.1  m/s . 
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The  smooth  domain  allows  a  variable  resolution  of 
grid  points  and  produces  excellent  results.  As  the  wave 
iensth  sets  shorter,  or  the  forecast  length  gets  larger, 
sere  s^all-scale  noise  is  apparent,  especially  with  the 
orimitive  f ormulat  ion  . 


3.   Arruot  Case 

The  abrupt  case  comparison  uses  the  abrupt  dorrain 
geometry  'see  Figure  16c).  This  grid  ooint  configuration  is 
v.sed  to  -"urther  illustrate  the  effects  of  spatial 
resolution.  The  previous  case  using  the  smooth  domain  and 
half  wave  length  introduced  noise  into  the  solution,  but  the 
spatial  resolution  changed  slowly  and  gradually. 

Consider  the  transition  necessary  in  an  operational 
model,  where  the  luxury  of  havi-n^-  a-  long  smooth-  transition 
into  the  region  of  high  resolution  may  not  be  possible.  The 
abrupt  domain  is  an  example  of  the  results  obtained  when 
spatial  resolution  is  decreased  rapidly. 

The  initial  fields  are  shown  in  Figure  26.  The  more 
linear  case,  A?A  =  1.1  m/s ,  is  used  to  eliminate  effects  due 
to  the  initial  field,  so  that  only  the  effects  due  to  the 
£riJ  ^ec^etry  are  seen. 

The  46  hour  comparisons  are  shown  in  Figure  27.  3oth 
solutions  are  affected  by  this  geometry,  but  the  primitive 
solutions  are  totally  disorganized  and  unacceptable. 

Table  4  shows  the  comparison  of  computational  times 
for  both  models  for  a  48  hour  forecast.  An   3££   savings   in 
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Figure  26.  Initial  fields  for  both  the  primitive  and  vorticity- 
divergence  models  using  the  abrupt  domain  and  perturbation  amplitude 
of  1.1  m/s.  Note  that  the  element  shapes  are  not  equilateral  triangles 
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Figure  27.  ^Q   hour  forecast  comparison  between  the  primitive  model 
and  the  vorticity-divergence  model  using  the  abrupt  domain.   (APA  = 
1.1  m/s). 
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CPU   tire   is  obtained  using  the  vortici ty-divergence  model. 
only  is   there  a   computational  savings   with   the 
vortici ty-divergence    model,    but    the   solution   is 
significantly  better  than  the  nrimitive  model. 


Table  4 

A   X 

'7v) 

c 
(m/s) 

At 
(sec) 

16? 
168 

320 

1? 

228 
6558 

vo^el    ax  c        At       #   of  steps     CPU  units 

■at  ions)      (sec) 

756  368.9 

25  55.7 

Comparison  of  the  computaional  times  for  a  48  hour  forecast 
between  the  primitive  model  and  the  vcrticity-di vergence 
roc" el  using  the  abrupt  domain. 


C.   COMPUTATIONAL  SENSITIVITY 

This  section  will  offer  an  explanation  for  the  skewing 
of  the  contours  in  the  96  h~our  forecast  solution- over  the 
smooth  domain  7i<?ure  23).  As  mentioned  previously,  there 
was  not  enough  time  available  for  fine  tuning  the 
c verrelaxation  coefficient,  which  is  used  while  solving  the 
system  of  equations.  The  overrelaxation  coefficient  may  be 
very  sensitive  and  small  changes  can,  on  occasion,  radically 
change  the  rate  of  convergence.  The  optimal  value  of  the 
overrelaxation  coefficient  depends  on  the  specific  form  of 
the  coefficient  of  the  equation  and  the  error  distribution. 

The  eouation  set  to  be  solved  consists  of  three 
equations  and  each  equation  required  its  own  relaxation 
coefficient.  When  solving  the  equations  over  the  regular 
domain,   the   entire   system   is  well  behaved  and  an  optimum 
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relaxation  coefficient  is  easily  determined.  However,  as  the 
ioirain  geometry  changes,  the   system  cf   equations   do   net 

converse  as  rapidly  and  the  relaxation  coefficient  requires 
further  fine  tuning. 

The  96  hour  forecast  uses  the  smooth  domain.  The 
f*i  d  -la  ti  tudinal  grid  points  are  compacted,  creates  a  denser 
belt  in  the  riddle  cf  the  channel.  The  coefficients 
originally  computed  using  the  regular  domain  need 
adjustments  to  properly  solve  the  equations. 

To  illustrate  the  significance  for  fine  tuning  the 
relaxation  coefficient,  consider  the  series  cf  plots  in 
Figure  28.  The  vorticity  divergence  equation  set  can  "be 
simplified  by  assuming  the  flow  is  non-divergent,  so  that 
only  the  vorticity  equation  needs  tc  be  solved.  Figure  23a 
is  the  46  hour  vorticity  field  using  this  equation  over  the 
regular  domain  with  an  o verrelaxa tion  coefficient  of  1.3. 
The  field  is  well  defined  with  smooth  contours. 

Figure  28b  is  similar  to  the  case  in  Figure  28a  except 
that  the  smooth  domain  is  used.  Notice  the  V-shaped  k:i nic  in 
the  pattern  wi*h  a  steeper  slope  in  the  upper  half.  This 
increased  bias  in  the  upper  half  is  caused  by  relaxing  the 
field  in  the  same  direction  during  each  pass  over  the 
domain.  When  the  direction  is  reversed  after  each  pass,  the 
exaggerated  bias  in  the  upper  half  disappears,  as  is  seen  in 
Figure  2£c.  However,  the  V-shaped  kink  is  net  eliminated. 
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Figure  28.  Computational  sensitivity  using  the  48  hour  vorticity  fields. 
Fig.  28  A)  the  48  hour  forecast  using  the  regular  domain,  28  B)  the  same 
forecast  using  the  smooth  domain  relaxing  in  one  direction,  28  C)  same  as 
28  B)  except  alternate  the  direction  of  relaxation  after  each  pass  over 
the  domain,  and  28  D)  same  as  28  C)  except  the  relaxation  coefficient  was 
fine  tuned  form  1.3  to  1.297. 
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Varying   the  relaxation   coefficient   from  1.3  to  1.297 
produces  the  much  improved  solution  in  Figure  28d  .   If   the 

relaxation  coefficients  for  the  other  two  equations  could 
alsa  he  fine  tuned,  it  appears  that  improved  solutions  would 
result.  There  are  also  other  relaxation  techniques  available 
that  have  potential  for  improving  the  solution  while  also 
sonrerging  at  a  faster  rate.  Some  of  these  techniques  will 
he  tested  in  the  futuer  usin^:  this  model. 
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VIII.  CONCLUSIONS 

T'his  research  investigated  different  finite  element 
formulations  for  the  shallow  water  equations.  The 
two— dimensional  dorrain  was  a  channel  which  simulated  a  belt 
around  the  earth.  Analytic  initial  conditions  were  used  to 
simplify  the  comparisons.  Two  formulations  were  examined; 
one  using  different  shaped  "basis  functions  and  the  other 
using  a  different  form  of  the  equation.  Zach  was  compared  to 
the  primitive  fornr  of  the  shallow  water  equations  that  was 
developed  by  Kelley  (1376). 

The  use  cf  equilateral  shaped  elements  which  was 
su^ested  by  ~r.  .V.J.F.  Cullen  significantly  improved  the 
solutions  compared  to  lelley's  model,  which  originally  used 
ri?ht  triangles  as  basis  functions,  *ost  of  the  other 
studies  in  this  thesis  used  the  equilateral  triangles. 

Williams  and  Zienkiewicz  (1981)  suggested  the  use  of 
piecewise  linear  basis  functions  for  the  velocity  field  and 
piecewise  constant  functions  for  the  height  field.  This 
formulation  was  tested  with  a  linearized  continuity 
equation.  The  results  were  poorer  than  those  obtained  with 
lelley's  model  . 

vost  of  the  effort  in  this  thesis  was  devoted  to 
ir-plement  at  ing   and   testing   a   vcrticity-divergence   model 


1  7C 
x  C  H 


sirrilar  to  the  ones  developed  by  Staniforth  and  Mitchell 
fig7?}  and  Cullen  and  Hall  (1979).  Several  tests  were 
presented  which  compare  this  formulation  with  Kelley's 
rod  el.  It  was  found  that  this  model  executes  approximately 
one  order  of  magnitude  faster  than  does  the  primitive 
formulation  ^sed  by  Kelley.  Secondly,  as  the  spatial 
resolution  between  grid  points  decreases,  this  formulation 
oroduces  a  solution  that  is  far  superior  to  the  primitive 
form.  A  disadvantage  is  its  computational  sensitivity,  which 
requires  fine  tuning  in  solving  the  elliptic  equations  for 
certain  geometries.  It  also  requires  25  percent  more 
computer  storage,  due  to  the  more  complex  equation  set  and 
the  additional  variables  that  are  treated. 

Implementation  of  finite  element  models  is  not  easy. 
Fowever,  there  is  a  lot  of  generality  and  redundancy 
imbedded  In  the  computations.  Versatile  modules  were  written 
which  significantly  reduced  the  effort  in  implementing  the 
vortici ty-d i vergence  model. 

further  research  is  suggested  using  this  finite  element 
formulation.  It  has  accurate  phase  propagation ,  is  able  to 
handle  variable  grid  geometry,  reduce  the  small-scale  noise 
and  decrease  the  model's  execution  time.  Specifically  more 
ad.va.r.cei  methods  of  solving  the  elliptic  equations  should  be 
investigated.  Finally,  the  formulation  should  be  tested  with 
small-scale  forcing,  where  its  advantages  should  be  most 
evident . 
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